diff --git a/algo/ca/core/data/CaWindowData.h b/algo/ca/core/data/CaWindowData.h index 0af05c83fcbde59432d0127632d4e4b16783927d..dc25c37a30f97c9afe275dae68979ee6198562cf 100644 --- a/algo/ca/core/data/CaWindowData.h +++ b/algo/ca/core/data/CaWindowData.h @@ -62,7 +62,7 @@ namespace cbm::algo::ca void ResetHitData(int nHits); /// \brief Reset suppressed hit flags - [[gnu::always_inline]] void ResetHitSuppressionFlags() { fvbHitSuppressed.reset(fvHits.size()); } + [[gnu::always_inline]] void ResetHitSuppressionFlags() { fvbHitSuppressed.reset(fvHits.size(), 0); } /// \brief Set hit suppression flag [[gnu::always_inline]] void SuppressHit(int iHit) { fvbHitSuppressed[iHit] = 1; } @@ -119,18 +119,18 @@ namespace cbm::algo::ca /// \brief Accesses reconstructed track container [[gnu::always_inline]] const Vector<Track>& RecoTracks() const { return fvRecoTracks; } - /// \brief Maps hit index from sub-TS to full dataset + /// \brief Maps hit index from the time window to the time slice /// \param iSt Index of tracking station - /// \param iHit Index of hit in the sub-TS - [[gnu::always_inline]] HitIndex_t FullDsHitIndex(int iSt, int iHit) const { return fvFullDsHitIndices[iSt][iHit]; } + /// \param iHit Index of hit in the window + [[gnu::always_inline]] HitIndex_t TsHitIndex(int iSt, int iHit) const { return fvTsHitIndices[iSt][iHit]; } - /// \brief Accesses container of hit index map from sub-TS to full dataset + /// \brief Accesses container of hit index map from the time window to the time slice /// \param iSt Index of tracking station - [[gnu::always_inline]] Vector<HitIndex_t>& FullDsHitIndices(int iSt) { return fvFullDsHitIndices[iSt]; } + [[gnu::always_inline]] Vector<HitIndex_t>& TsHitIndices(int iSt) { return fvTsHitIndices[iSt]; } - /// \brief Accesses container of hit index map from sub-TS to full dataset + /// \brief Accesses container of hit index map from the time window to the time slice /// \param iSt Index of tracking station - [[gnu::always_inline]] const Vector<HitIndex_t>& FullDsHitIndices(int iSt) const { return fvFullDsHitIndices[iSt]; } + [[gnu::always_inline]] const Vector<HitIndex_t>& TsHitIndices(int iSt) const { return fvTsHitIndices[iSt]; } /// \brief Accesses current iteration [[gnu::always_inline]] void SetCurrentIteration(const Iteration* ptr) { fpCurrentIteration = ptr; } @@ -138,19 +138,17 @@ namespace cbm::algo::ca /// \brief Accesses current iteration [[gnu::always_inline]] const Iteration* CurrentIteration() const { return fpCurrentIteration; } - fvec& TargX() { return fTargX; } - fvec& TargY() { return fTargY; } - fvec& TargZ() { return fTargZ; } + /// \brief Accesses magnetic field in starting point (target or first station) + [[gnu::always_inline]] ca::FieldValue<fvec>& TargB() { return fTargB; } - fvec TargX() const { return fTargX; } - fvec TargY() const { return fTargY; } - fvec TargZ() const { return fTargZ; } + /// \brief Accesses magnetic field in starting point (target or first station) + [[gnu::always_inline]] const ca::FieldValue<fvec>& TargB() const { return fTargB; } - ca::FieldValue<fvec>& TargB() { return fTargB; } - const ca::FieldValue<fvec>& TargB() const { return fTargB; } + /// \brief Measurement of the target with the uncertainty + [[gnu::always_inline]] ca::MeasurementXy<fvec>& TargetMeasurement() { return fTargetMeasurement; } - ca::MeasurementXy<fvec>& TargetMeasurement() { return fTargetMeasurement; } - const ca::MeasurementXy<fvec>& TargetMeasurement() const { return fTargetMeasurement; } + /// \brief Measurement of the target with the uncertainty + [[gnu::always_inline]] const ca::MeasurementXy<fvec>& TargetMeasurement() const { return fTargetMeasurement; } private: static constexpr int kMaxNofStations = constants::size::MaxNstations; ///< Alias to max number of stations @@ -182,18 +180,14 @@ namespace cbm::algo::ca /// \brief Sample of track candidates Vector<Track> fvTrackCandidates{"WindowData::fvTrackCandidates"}; - /// \brief Map of hit indices from sub-TS to full dataset - std::array<Vector<HitIndex_t>, kMaxNofStations> fvFullDsHitIndices{"WindowData::fvFullDSHitIndex"}; + /// \brief Map of hit indices from the time window to the time slice + std::array<Vector<HitIndex_t>, kMaxNofStations> fvTsHitIndices{"WindowData::fvFullDSHitIndex"}; /// \brief Current track-finder iteration const Iteration* fpCurrentIteration = nullptr; - fvec fTargX{constants::Undef<fvec>}; ///< target position x coordinate for the current iteration (modifiable) - fvec fTargY{constants::Undef<fvec>}; ///< target position y coordinate for the current iteration (modifiable) - fvec fTargZ{constants::Undef<fvec>}; ///< target position z coordinate for the current iteration (modifiable) - - ca::FieldValue<fvec> fTargB _fvecalignment{}; // field in the target point (modifiable, do not touch!!) - ca::MeasurementXy<fvec> fTargetMeasurement _fvecalignment{}; // target constraint [cm] + ca::FieldValue<fvec> fTargB _fvecalignment{}; ///< field in the target point (modifiable, do not touch!!) + ca::MeasurementXy<fvec> fTargetMeasurement _fvecalignment{}; ///< target constraint } _fvecalignment; } // namespace cbm::algo::ca diff --git a/algo/ca/core/pars/CaIteration.h b/algo/ca/core/pars/CaIteration.h index 0704d3c103950f4e9d21bb79e0fc688d2de71869..1ec0e78a5c29ab666f64d7b86f9c7764cb219218 100644 --- a/algo/ca/core/pars/CaIteration.h +++ b/algo/ca/core/pars/CaIteration.h @@ -144,12 +144,7 @@ namespace cbm::algo::ca /// \brief Sets flag: true - triplets are built also skipping <= GetMaxStationGap stations void SetMaxStationGap(int nSkipped) { fMaxStationGap = nSkipped; } - /// \brief Sets correction for accounting overlaping and iff z - /// - /// the maxDZ correction is set in order to take into account overlaping and iff z. - /// The reason is that low momentum tracks are too curved and goes not from target direction. - /// That's why sort by hit_y/hit_z is not work idealy. If sort by y then it is max diff between - /// same station's modules (~0.4cm). + /// \brief TODO: select a more proper name void SetMaxDZ(float input) { fMaxDZ = input; } /// \brief Sets max considered q/p for tracks diff --git a/algo/ca/core/tracking/CaCloneMerger.h b/algo/ca/core/tracking/CaCloneMerger.h index 38e4d89def5fee76166f1f744c85a4a075dc5c40..c16d53fe4fed6220d2552f287f0633d05128f31c 100644 --- a/algo/ca/core/tracking/CaCloneMerger.h +++ b/algo/ca/core/tracking/CaCloneMerger.h @@ -36,10 +36,10 @@ namespace cbm::algo::ca ~CloneMerger(); /// Copy constructor - CloneMerger(const CloneMerger&) = delete; + CloneMerger(const CloneMerger&) = default; /// Move constructor - CloneMerger(CloneMerger&&) = delete; + CloneMerger(CloneMerger&&) = default; /// Copy assignment operator CloneMerger& operator=(const CloneMerger&) = delete; diff --git a/algo/ca/core/tracking/CaFramework.cxx b/algo/ca/core/tracking/CaFramework.cxx index 96907d4dec62c9d0395733af21ecfd8dbe2d02f4..cc3ebd1121507493ea532a84d47d9e2b54af0d4a 100644 --- a/algo/ca/core/tracking/CaFramework.cxx +++ b/algo/ca/core/tracking/CaFramework.cxx @@ -24,8 +24,23 @@ using cbm::algo::ca::EDetectorID; using cbm::algo::ca::ETimer; // monitor timer key type //using cbm::ca::tools::Debugger; -void Framework::Init(const TrackingMode mode) { fTrackingMode = mode; } +// --------------------------------------------------------------------------------------------------------------------- +// +void Framework::Init(const TrackingMode mode) +{ + fTrackingMode = mode; + for (int iThread = 0; iThread < fNofThreads; ++iThread) { + fvWData.emplace_back(); + LOG(info) << "->>>" << &fvWData.back(); + } + for (int iThread = 0; iThread < fNofThreads; ++iThread) { + fvTrackFinderWindow.emplace_back(*this, fvWData[iThread]); + } + fTrackFinder.Init(); +} +// --------------------------------------------------------------------------------------------------------------------- +// void Framework::Finish() { //Debugger::Instance().Write(); diff --git a/algo/ca/core/tracking/CaFramework.h b/algo/ca/core/tracking/CaFramework.h index 0d4c18ad3c2c7d663dcd41f37cf55ac677b70e5f..50e1a7069f01a47f4f1f27da8cc674198e2e93ec 100644 --- a/algo/ca/core/tracking/CaFramework.h +++ b/algo/ca/core/tracking/CaFramework.h @@ -199,6 +199,15 @@ namespace cbm::algo::ca // const CbmL1MCTrack* GetMcTrackForWindowHit(int iHit) const; + void SetNofThtreads(int nThreads) + { + fNofThreads = nThreads; + LOG(info) << "ca::Framework: number of threads is set to " << fNofThreads; + } + + int GetNofThreads() const { return fNofThreads; } + + private: int fNstationsBeforePipe{0}; ///< number of stations before pipe (MVD stations in CBM) int fNfieldStations{0}; ///< number of stations in the field region @@ -217,6 +226,8 @@ namespace cbm::algo::ca TrackingMonitorData fMonitorData{}; ///< Tracking monitor data (statistics per call) + int fNofThreads = 1; ///< Number of threads to execute the track-finder + public: Vector<CaHitTimeInfo> fHitTimeInfo; @@ -236,7 +247,6 @@ namespace cbm::algo::ca // WARN: Potential race conditions -> Vector<int> fHitKeyToTrack{"Framework::fHitKeyToTrack"}; // strip to track pointers - TrackingMode fTrackingMode{kSts}; fvec EventTime{0.f}; fvec Err{0.f}; @@ -245,12 +255,10 @@ namespace cbm::algo::ca unsigned int fCurrentIterationIndex{0}; ///< index of the corrent iteration, needed for debug output public: - ca::WindowData fWData{}; - ca::TrackExtender fTrackExtender{*this}; ///< Object of the track extender algorithm - ca::CloneMerger fCloneMerger{*this}; ///< Object of the clone merger algorithm - ca::TrackFitter fTrackFitter{*this}; ///< Object of the track extender algorithm - ca::TrackFinderWindow fTrackFinderWindow{*this}; ///< Object of the track finder algorithm for the time window - ca::TrackFinder fTrackFinder{*this}; ///< Object of the track finder algorithm for the entire time slice + 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 ================================= diff --git a/algo/ca/core/tracking/CaTrackExtender.cxx b/algo/ca/core/tracking/CaTrackExtender.cxx index efd6418ea3f3e6398331b11907eb183654ac4918..9dcfddbd2a4a2e9e2559b7aae586dd2201eec5d6 100644 --- a/algo/ca/core/tracking/CaTrackExtender.cxx +++ b/algo/ca/core/tracking/CaTrackExtender.cxx @@ -23,7 +23,7 @@ using namespace cbm::algo::ca; // --------------------------------------------------------------------------------------------------------------------- // -TrackExtender::TrackExtender(const ca::Framework& algo) : frAlgo(algo) {} +TrackExtender::TrackExtender(const ca::Framework& algo, WindowData& wData) : frAlgo(algo), frWData(wData) {} // --------------------------------------------------------------------------------------------------------------------- // @@ -192,11 +192,11 @@ 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 == frAlgo.fWData.CurrentIteration()->GetFirstStationIndex()) ista = ista2 + step; + if (ista2 == frWData.CurrentIteration()->GetFirstStationIndex()) ista = ista2 + step; - const float pickGather = frAlgo.fWData.CurrentIteration()->GetPickGather(); + const float pickGather = frWData.CurrentIteration()->GetPickGather(); const fvec pickGather2 = pickGather * pickGather; - const fvec maxDZ = frAlgo.fWData.CurrentIteration()->GetMaxDZ(); + const fvec maxDZ = frWData.CurrentIteration()->GetMaxDZ(); for (; (ista < frAlgo.GetParameters().GetNstationsActive()) && (ista >= 0); ista += step) { // CHECKME why ista2? const ca::Station<fvec>& sta = frAlgo.GetParameters().GetStation(ista); @@ -208,7 +208,7 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up TrackParamV& tr = fit.Tr(); - const auto& grid = frAlgo.fWData.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]); @@ -221,10 +221,10 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up while (area.GetNextObjectId(ih)) { // loop over the hits in the area - if (frAlgo.fWData.IsHitSuppressed(ih)) { + if (frWData.IsHitSuppressed(ih)) { continue; } - const ca::Hit& hit = frAlgo.fWData.Hit(ih); + const ca::Hit& hit = frWData.Hit(ih); if (sta.timeInfo && tr.NdfTime()[0] > -2.) { fscal dt = hit.T() - tr.Time()[0]; @@ -262,7 +262,7 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up if (iHit_best < 0) break; - const ca::Hit& hit = frAlgo.fWData.Hit(iHit_best); + const ca::Hit& hit = frWData.Hit(iHit_best); newHits.push_back(hit.Id()); fit.Extrapolate(hit.Z(), fld); diff --git a/algo/ca/core/tracking/CaTrackExtender.h b/algo/ca/core/tracking/CaTrackExtender.h index fe14ed39c6eeb4539fb0114302fce55290e799e1..e5762fb67e67214d32146c209a77d15c389e7c97 100644 --- a/algo/ca/core/tracking/CaTrackExtender.h +++ b/algo/ca/core/tracking/CaTrackExtender.h @@ -14,6 +14,7 @@ #include "CaSimd.h" #include "CaTrackParam.h" #include "CaVector.h" +#include "CaWindowData.h" namespace cbm::algo::ca @@ -32,16 +33,16 @@ namespace cbm::algo::ca class TrackExtender { public: /// Default constructor - TrackExtender(const ca::Framework& algo); + TrackExtender(const ca::Framework& algo, WindowData& wData); /// Destructor ~TrackExtender(); /// Copy constructor - TrackExtender(const TrackExtender&) = delete; + TrackExtender(const TrackExtender&) = default; /// Move constructor - TrackExtender(TrackExtender&&) = delete; + TrackExtender(TrackExtender&&) = default; /// Copy assignment operator TrackExtender& operator=(const TrackExtender&) = delete; @@ -88,6 +89,7 @@ namespace cbm::algo::ca /// Data members const ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class + WindowData& frWData; }; } // namespace cbm::algo::ca diff --git a/algo/ca/core/tracking/CaTrackFinder.cxx b/algo/ca/core/tracking/CaTrackFinder.cxx index c466fa1fd82bec47ec4af675afc83eeb9b7482b7..c997f2fa1a669d69f2e1d5a0efdf87e993c29a5e 100644 --- a/algo/ca/core/tracking/CaTrackFinder.cxx +++ b/algo/ca/core/tracking/CaTrackFinder.cxx @@ -61,9 +61,11 @@ void TrackFinder::FindTracks() frAlgo.fRecoHits.reserve(2 * frAlgo.fInputData.GetNhits()); frAlgo.fRecoTracks.reserve(2 * frAlgo.fInputData.GetNhits() / frAlgo.GetParameters().GetNstationsActive()); - for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); ++iS) { - frAlgo.fWData.FullDsHitIndices(iS).clear(); - frAlgo.fWData.FullDsHitIndices(iS).reserve(frAlgo.fInputData.GetNhits()); + for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) { + for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); ++iS) { + frAlgo.fvWData[iThread].TsHitIndices(iS).clear(); + frAlgo.fvWData[iThread].TsHitIndices(iS).reserve(frAlgo.fInputData.GetNhits()); + } } frAlgo.fvHitKeyFlags.reset(frAlgo.fInputData.GetNhitKeys(), 0); @@ -74,9 +76,6 @@ void TrackFinder::FindTracks() // TODO: move these values to Parameters namespace (S.Zharko) // length of sub-TS - const fscal tsLength = (ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) ? 500 : 10000; - - const fscal tsOverlap = 15; // min length of overlap region fscal minProtonMomentum = 0.1; fscal targX = frAlgo.GetParameters().GetTargetPositionX()[0]; @@ -131,7 +130,7 @@ void TrackFinder::FindTracks() if (info.fEventTimeMin > 500.e6 || info.fEventTimeMax < -500.) { // cut hits with bogus start time > 500 ms frAlgo.fvHitKeyFlags[h.FrontKey()] = 1; frAlgo.fvHitKeyFlags[h.BackKey()] = 1; - LOG(warning) << "CA: skip bogus hit " << h.ToString(); + LOG(error) << "CATrackFinder: skip bogus hit " << h.ToString(); continue; } @@ -202,153 +201,151 @@ void TrackFinder::FindTracks() // cut data into sub-timeslices and process them one by one - bool areUntouchedDataLeft = true; // is the whole TS processed + for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) { + bool areUntouchedDataLeft = true; // is the whole TS processed + ca::HitIndex_t sliceFirstHit[nDataStreams]; - ca::HitIndex_t sliceFirstHit[nDataStreams]; + for (int iStream = 0; iStream < nDataStreams; ++iStream) { + sliceFirstHit[iStream] = frAlgo.fInputData.GetStreamStartIndex(iStream); + } - for (int iStream = 0; iStream < nDataStreams; ++iStream) { - sliceFirstHit[iStream] = frAlgo.fInputData.GetStreamStartIndex(iStream); - } + int statNwindows = 0; + int statNhitsProcessed = 0; + int statLastLogTimeChunk = -1; - int statNwindows = 0; - int statNhitsProcessed = 0; - int statLastLogTimeChunk = -1; + while (areUntouchedDataLeft) { + frAlgo.fMonitorData.IncrementCounter(ECounter::SubTS); + // select the sub-slice hits + for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); ++iS) { + frAlgo.fvWData[iThread].TsHitIndices(iS).clear(); + } - while (areUntouchedDataLeft) { - frAlgo.fMonitorData.IncrementCounter(ECounter::SubTS); - // select the sub-slice hits - for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); ++iS) { - frAlgo.fWData.FullDsHitIndices(iS).clear(); - } + areUntouchedDataLeft = false; - areUntouchedDataLeft = false; + // TODO: SG: skip empty regions and start the subslice with the earliest hit - // TODO: SG: skip empty regions and start the subslice with the earliest hit + statNwindows++; + int statNwindowHits = 0; - statNwindows++; - int statNwindowHits = 0; + for (int iStream = 0; iStream < nDataStreams; ++iStream) { - for (int iStream = 0; iStream < nDataStreams; ++iStream) { + for (ca::HitIndex_t caHitId = sliceFirstHit[iStream]; caHitId < frAlgo.fInputData.GetStreamStopIndex(iStream); + ++caHitId) { - for (ca::HitIndex_t caHitId = sliceFirstHit[iStream]; caHitId < frAlgo.fInputData.GetStreamStopIndex(iStream); - ++caHitId) { - - CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId]; - const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId); + CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId]; + const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId); - if (frAlgo.fvHitKeyFlags[h.FrontKey()] - || frAlgo.fvHitKeyFlags[h.BackKey()]) { // the hit is already reconstructed - continue; - } - if (info.fEventTimeMin > tsStart + tsLength) { // the hit is too late for the sub slice - areUntouchedDataLeft = true; - if (info.fMinTimeAfterHit > tsStart + tsLength) { - // this hit and all later hits are out of the sub-slice - break; + if (frAlgo.fvHitKeyFlags[h.FrontKey()] + || frAlgo.fvHitKeyFlags[h.BackKey()]) { // the hit is already reconstructed + continue; } - } - else { - if (tsStart <= info.fEventTimeMax) { // the hit belongs to the sub-slice - frAlgo.fWData.FullDsHitIndices(h.Station()).push_back(caHitId); - statNwindowHits++; - if (info.fMaxTimeBeforeHit < tsStart + tsLength - tsOverlap) { - // this hit and all hits before are before the overlap - sliceFirstHit[iStream] = caHitId + 1; + if (info.fEventTimeMin > tsStart + fWindowLength) { // the hit is too late for the sub slice + areUntouchedDataLeft = true; + if (info.fMinTimeAfterHit > tsStart + fWindowLength) { + // this hit and all later hits are out of the sub-slice + break; } - //LOG(warning) << "event time " << ((float) info.fEventTimeMax )<< " " << h.ToString(); } - } // else the hit has been alread processed in previous sub-slices + else { + if (tsStart <= info.fEventTimeMax) { // the hit belongs to the sub-slice + frAlgo.fvWData[iThread].TsHitIndices(h.Station()).push_back(caHitId); + statNwindowHits++; + if (info.fMaxTimeBeforeHit < tsStart + fWindowLength - fWindowOverlap) { + // this hit and all hits before are before the overlap + sliceFirstHit[iStream] = caHitId + 1; + } + //LOG(warning) << "event time " << ((float) info.fEventTimeMax )<< " " << h.ToString(); + } + } // else the hit has been alread processed in previous sub-slices + } } - } - statNhitsProcessed += statNwindowHits; - - // print the LOG for every 10 ms of data processed - { - int currentChunk = (int) ((tsStart - statTsStart) / 10.e6); - if (!areUntouchedDataLeft || currentChunk > statLastLogTimeChunk) { - statLastLogTimeChunk = currentChunk; - double dataRead = 100. * (tsStart + tsLength - statTsStart) / (statTsEnd - statTsStart); - if (dataRead > 100.) { - dataRead = 100.; + statNhitsProcessed += statNwindowHits; + + // print the LOG for every 10 ms of data processed + { + int currentChunk = (int) ((tsStart - statTsStart) / 10.e6); + if (!areUntouchedDataLeft || currentChunk > statLastLogTimeChunk) { + statLastLogTimeChunk = currentChunk; + double dataRead = 100. * (tsStart + fWindowLength - statTsStart) / (statTsEnd - statTsStart); + if (dataRead > 100.) { + dataRead = 100.; + } + LOG(info) << "CA tracker process sliding window N " << statNwindows << ": time " << tsStart / 1.e6 << " ms + " + << fWindowLength / 1.e3 << " us) with " << statNwindowHits << " hits. " + << " Processing " << dataRead << " % of the TS time and " + << 100. * statNhitsProcessed / statNhitsTotal << " % of TS hits." + << " Already reconstructed " << frAlgo.fRecoTracks.size() << " tracks "; } - LOG(info) << "CA tracker process sliding window N " << statNwindows << ": time " << tsStart / 1.e6 << " ms + " - << tsLength / 1.e3 << " us) with " << statNwindowHits << " hits. " - << " Processing " << dataRead << " % of the TS time and " - << 100. * statNhitsProcessed / statNhitsTotal << " % of TS hits." - << " Already reconstructed " << frAlgo.fRecoTracks.size() << " tracks "; } - } - if (ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) { - // cut at 50 hits per station per 1 us. - int maxStationHits = (int) (50 * tsLength / 1.e3); - for (int ista = 0; ista < frAlgo.GetParameters().GetNstationsActive(); ++ista) { - if ((int) frAlgo.fWData.FullDsHitIndices(ista).size() > maxStationHits) { - frAlgo.fWData.FullDsHitIndices(ista).clear(); + if (ca::Framework::TrackingMode::kMcbm == frAlgo.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) { + if ((int) frAlgo.fvWData[iThread].TsHitIndices(ista).size() > maxStationHits) { + frAlgo.fvWData[iThread].TsHitIndices(ista).clear(); + } } } - } - frAlgo.fMonitorData.StartTimer(ETimer::TrackFinder); - frAlgo.fTrackFinderWindow.CaTrackFinderSlice(); - frAlgo.fMonitorData.StopTimer(ETimer::TrackFinder); + frAlgo.fMonitorData.StartTimer(ETimer::TrackFinder); + frAlgo.fvTrackFinderWindow[iThread].CaTrackFinderSlice(); + frAlgo.fMonitorData.StopTimer(ETimer::TrackFinder); - frAlgo.fMonitorData.StartTimer(ETimer::TrackFitter); - frAlgo.fTrackFitter.FitCaTracks(); - frAlgo.fMonitorData.StopTimer(ETimer::TrackFitter); - // save reconstructed tracks with no hits in the overlap region + // save reconstructed tracks with no hits in the overlap region - tsStart += tsLength - tsOverlap; + tsStart += 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 + // 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 - int trackFirstHit = 0; - for (const auto& track : frAlgo.fWData.RecoTracks()) { - bool isTrackCompletelyInOverlap = true; - for (int ih = 0; ih < track.fNofHits; ih++) { - int caHitId = frAlgo.fWData.RecoHitIndex(trackFirstHit + ih); - CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId]; - if (info.fEventTimeMax < tsStart) { // this hit is before the overlap - isTrackCompletelyInOverlap = false; - break; + 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 < tsStart) { // this hit is before the overlap + isTrackCompletelyInOverlap = false; + break; + } } - } - if (!areUntouchedDataLeft) { // don't reject tracks in the overlap when no more data are left - isTrackCompletelyInOverlap = 0; - } + if (!areUntouchedDataLeft) { // don't reject tracks in the overlap when no more data are left + isTrackCompletelyInOverlap = 0; + } - 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.fWData.RecoHitIndex(trackFirstHit + i); - const auto& h = frAlgo.fInputData.GetHit(caHitId); - frAlgo.fvHitKeyFlags[h.FrontKey()] = 0; - frAlgo.fvHitKeyFlags[h.BackKey()] = 0; + 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); + frAlgo.fvHitKeyFlags[h.FrontKey()] = 0; + frAlgo.fvHitKeyFlags[h.BackKey()] = 0; + } } - } - else { // save the track - frAlgo.fRecoTracks.push_back(track); - // mark the track hits as used - for (int i = 0; i < track.fNofHits; i++) { - int caHitId = frAlgo.fWData.RecoHitIndex(trackFirstHit + i); - const auto& h = frAlgo.fInputData.GetHit(caHitId); - frAlgo.fvHitKeyFlags[h.FrontKey()] = 1; - frAlgo.fvHitKeyFlags[h.BackKey()] = 1; - frAlgo.fRecoHits.push_back(caHitId); + else { // save the track + frAlgo.fRecoTracks.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); + frAlgo.fvHitKeyFlags[h.FrontKey()] = 1; + frAlgo.fvHitKeyFlags[h.BackKey()] = 1; + frAlgo.fRecoHits.push_back(caHitId); + } } - } - trackFirstHit += track.fNofHits; - } // sub-timeslice tracks + trackFirstHit += track.fNofHits; + } // sub-timeslice tracks - tsStart -= 5; // do 5 ns margin + tsStart -= 5; // do 5 ns margin + } } frAlgo.fMonitorData.StopTimer(ETimer::FindTracks); @@ -358,3 +355,18 @@ void TrackFinder::FindTracks() auto timerEnd = std::chrono::high_resolution_clock::now(); frAlgo.fCaRecoTime = (double) (std::chrono::duration<double>(timerEnd - timerStart).count()); } + +// --------------------------------------------------------------------------------------------------------------------- +// +void TrackFinder::Init() +{ + fvTsStartThread.clear(); + fvTsEndThread.clear(); + fvWindowStartThread.clear(); + + fvTsStartThread.resize(frAlgo.GetNofThreads()); + fvTsEndThread.resize(frAlgo.GetNofThreads()); + fvWindowStartThread.resize(frAlgo.GetNofThreads()); + + 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 71420a371379ddb2b0d482402c106a58eed069f2..9bfe0f9bcb77cc3c8f3b935d1bb9e66d92c20cef 100644 --- a/algo/ca/core/tracking/CaTrackFinder.h +++ b/algo/ca/core/tracking/CaTrackFinder.h @@ -48,6 +48,7 @@ namespace cbm::algo::ca TrackFinder& operator=(TrackFinder&&) = delete; void FindTracks(); + void Init(); private: // ------------------------------- @@ -61,6 +62,13 @@ namespace cbm::algo::ca // Data members ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class + + std::vector<float> fvTsStartThread; + std::vector<float> fvTsEndThread; + std::vector<float> fvWindowStartThread; + + float fWindowLength = 0.; + float fWindowOverlap = 15.; // ns }; } // namespace cbm::algo::ca diff --git a/algo/ca/core/tracking/CaTrackFinderWindow.cxx b/algo/ca/core/tracking/CaTrackFinderWindow.cxx index 0336d2b10c0489bd7eb64511d9d92003221cd7af..f25135dc4c5ea703a04e6299231947a67c3cefb1 100644 --- a/algo/ca/core/tracking/CaTrackFinderWindow.cxx +++ b/algo/ca/core/tracking/CaTrackFinderWindow.cxx @@ -46,13 +46,19 @@ using namespace cbm::algo::ca; // --------------------------------------------------------------------------------------------------------------------- // -TrackFinderWindow::TrackFinderWindow(ca::Framework& algo) : frAlgo(algo) +TrackFinderWindow::TrackFinderWindow(ca::Framework& algo, ca::WindowData& wData) + : frAlgo(algo) + , frWData(wData) + , fTrackExtender(frAlgo, frWData) + , fCloneMerger(frAlgo) + , fTrackFitter(frAlgo, frWData) { for (unsigned int i = 0; i < constants::size::MaxNstations; ++i) { fvTriplets[i].SetName(std::stringstream() << "TrackFinderWindow::fTriplets[" << i << "]"); } } + // --------------------------------------------------------------------------------------------------------------------- // TrackFinderWindow::~TrackFinderWindow() {} @@ -72,7 +78,7 @@ bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triple if (r.GetLSta() != l.GetMSta()) return false; - const float tripletLinkChi2 = frAlgo.fWData.CurrentIteration()->GetTripletLinkChi2(); + const float tripletLinkChi2 = frWData.CurrentIteration()->GetTripletLinkChi2(); if (r.IsMomentumFitted()) { assert(l.IsMomentumFitted()); @@ -126,26 +132,26 @@ void TrackFinderWindow::ReadWindowData() { int nHits = 0; for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); iS++) { - frAlgo.fWData.HitStartIndexOnStation(iS) = nHits; - frAlgo.fWData.NofHitsOnStation(iS) = frAlgo.fWData.FullDsHitIndices(iS).size(); - nHits += frAlgo.fWData.NofHitsOnStation(iS); + frWData.HitStartIndexOnStation(iS) = nHits; + frWData.NofHitsOnStation(iS) = frWData.TsHitIndices(iS).size(); + nHits += frWData.NofHitsOnStation(iS); } - frAlgo.fWData.HitStartIndexOnStation(frAlgo.GetParameters().GetNstationsActive()) = nHits; - frAlgo.fWData.ResetHitData(nHits); + frWData.HitStartIndexOnStation(frAlgo.GetParameters().GetNstationsActive()) = nHits; + frWData.ResetHitData(nHits); for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); iS++) { - int iFstHit = frAlgo.fWData.HitStartIndexOnStation(iS); - for (ca::HitIndex_t ih = 0; ih < frAlgo.fWData.FullDsHitIndices(iS).size(); ++ih) { - ca::Hit h = frAlgo.fInputData.GetHit(frAlgo.fWData.FullDsHitIndex(iS, ih)); - h.SetId(frAlgo.fWData.FullDsHitIndex(iS, ih)); - frAlgo.fWData.Hit(iFstHit + ih) = h; + 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; } } if constexpr (fDebug) { LOG(info) << "===== Sliding Window hits: "; for (int i = 0; i < nHits; ++i) { - LOG(info) << " " << frAlgo.fWData.Hit(i).ToString(); + LOG(info) << " " << frWData.Hit(i).ToString(); } LOG(info) << "===== "; } @@ -154,16 +160,16 @@ void TrackFinderWindow::ReadWindowData() fvHitNofTriplets.reset(nHits, 0); - frAlgo.fWData.RecoTracks().clear(); - frAlgo.fWData.RecoTracks().reserve(2 * nHits / frAlgo.GetParameters().GetNstationsActive()); + frWData.RecoTracks().clear(); + frWData.RecoTracks().reserve(2 * nHits / frAlgo.GetParameters().GetNstationsActive()); - frAlgo.fWData.RecoHitIndices().clear(); - frAlgo.fWData.RecoHitIndices().reserve(2 * nHits); + frWData.RecoHitIndices().clear(); + frWData.RecoHitIndices().reserve(2 * nHits); fvTrackCandidates.clear(); fvTrackCandidates.reserve(nHits / 10); for (unsigned int iS = 0; iS < constants::size::MaxNstations; iS++) { - int nHitsStation = frAlgo.fWData.FullDsHitIndices(iS).size(); + int nHitsStation = frWData.TsHitIndices(iS).size(); fvTriplets[iS].clear(); fvTriplets[iS].reserve(2 * nHitsStation); } @@ -189,8 +195,8 @@ void TrackFinderWindow::CaTrackFinderSlice() fscal gridMinY = -0.1; fscal gridMaxY = 0.1; - for (ca::HitIndex_t ih = 0; ih < frAlgo.fWData.FullDsHitIndices(iS).size(); ++ih) { - const ca::Hit& h = frAlgo.fInputData.GetHit(frAlgo.fWData.FullDsHitIndex(iS, ih)); + for (ca::HitIndex_t ih = 0; ih < frWData.TsHitIndices(iS).size(); ++ih) { + const ca::Hit& h = frAlgo.fInputData.GetHit(frWData.TsHitIndex(iS, ih)); if (h.X() < gridMinX) { gridMinX = h.X(); @@ -219,7 +225,7 @@ void TrackFinderWindow::CaTrackFinderSlice() // TODO: changing the grid also changes the result. Investigate why it happens. // TODO: find the optimal grid size - int nSliceHits = frAlgo.fWData.FullDsHitIndices(iS).size(); + int nSliceHits = frWData.TsHitIndices(iS).size(); fscal sizeY = gridMaxY - gridMinY; fscal sizeX = gridMaxX - gridMinX; int nBins2D = 1 + nSliceHits; @@ -234,12 +240,12 @@ void TrackFinderWindow::CaTrackFinderSlice() if (xStep > 0.3 * scale) xStep = 0.3 * scale; if (xStep < 0.01 * scale) xStep = 0.01 * scale; - auto& grid = frAlgo.fWData.Grid(iS); + auto& grid = frWData.Grid(iS); grid.BuildBins(gridMinX, gridMaxX, gridMinY, gridMaxY, xStep, yStep); /* clang-format off */ - grid.StoreHits(frAlgo.fWData.Hits(), - frAlgo.fWData.HitStartIndexOnStation(iS), - frAlgo.fWData.NofHitsOnStation(iS), + grid.StoreHits(frWData.Hits(), + frWData.HitStartIndexOnStation(iS), + frWData.NofHitsOnStation(iS), frAlgo.fvHitKeyFlags); /* clang-format on */ } @@ -253,17 +259,17 @@ void TrackFinderWindow::CaTrackFinderSlice() // current iteration of the CA tracker const auto& caIteration = frAlgo.GetParameters().GetCAIterations()[frAlgo.fCurrentIterationIndex]; - frAlgo.fWData.SetCurrentIteration(&caIteration); + frWData.SetCurrentIteration(&caIteration); if (frAlgo.fCurrentIterationIndex > 0) { for (int ista = 0; ista < frAlgo.GetParameters().GetNstationsActive(); ++ista) { - frAlgo.fWData.Grid(ista).RemoveUsedHits(frAlgo.fWData.Hits(), frAlgo.fvHitKeyFlags); + frWData.Grid(ista).RemoveUsedHits(frWData.Hits(), frAlgo.fvHitKeyFlags); } } - int nHits = frAlgo.fWData.Hits().size(); + int nHits = frWData.Hits().size(); // --> frAlgo.fIsWindowHitSuppressed.reset(frAlgo.fWindowHits.size(), 0); - frAlgo.fWData.ResetHitSuppressionFlags(); // TODO: ??? No effect? + frWData.ResetHitSuppressionFlags(); // TODO: ??? No effect? for (int j = 0; j < frAlgo.GetParameters().GetNstationsActive(); j++) { fvTriplets[j].clear(); } @@ -276,30 +282,26 @@ void TrackFinderWindow::CaTrackFinderSlice() { // --- SET PARAMETERS FOR THE ITERATION --- // define the target - frAlgo.fWData.TargX() = frAlgo.GetParameters().GetTargetPositionX(); - frAlgo.fWData.TargY() = frAlgo.GetParameters().GetTargetPositionY(); - frAlgo.fWData.TargZ() = frAlgo.GetParameters().GetTargetPositionZ(); - fscal SigmaTargetX = caIteration.GetTargetPosSigmaX(); fscal SigmaTargetY = caIteration.GetTargetPosSigmaY(); // target constraint [cm] // Select magnetic field. For primary tracks - fVtxFieldValue, for secondary tracks - st.fieldSlice if (caIteration.GetPrimaryFlag()) { - frAlgo.fWData.TargB() = frAlgo.GetParameters().GetVertexFieldValue(); + frWData.TargB() = frAlgo.GetParameters().GetVertexFieldValue(); } else { - frAlgo.fWData.TargB() = frAlgo.GetParameters().GetStation(0).fieldSlice.GetFieldValue(0, 0); + frWData.TargB() = frAlgo.GetParameters().GetStation(0).fieldSlice.GetFieldValue(0, 0); } // NOTE: calculates field frAlgo.fTargB in the center of 0th station - frAlgo.fIsTargetField = (fabs(frAlgo.fWData.TargB().y[0]) > 0.001); + frAlgo.fIsTargetField = (fabs(frWData.TargB().y[0]) > 0.001); - frAlgo.fWData.TargetMeasurement().SetX(frAlgo.fWData.TargX()); - frAlgo.fWData.TargetMeasurement().SetY(frAlgo.fWData.TargY()); - frAlgo.fWData.TargetMeasurement().SetDx2(SigmaTargetX * SigmaTargetX); - frAlgo.fWData.TargetMeasurement().SetDxy(0); - frAlgo.fWData.TargetMeasurement().SetDy2(SigmaTargetY * SigmaTargetY); - frAlgo.fWData.TargetMeasurement().SetNdfX(1); - frAlgo.fWData.TargetMeasurement().SetNdfY(1); + 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); } } @@ -308,20 +310,20 @@ void TrackFinderWindow::CaTrackFinderSlice() frAlgo.fMonitorData.StartTimer(ETimer::TripletConstruction); - ca::TripletConstructor constructor(&frAlgo); + ca::TripletConstructor constructor(frAlgo, frWData); // indices of the two neighbouring station, taking into account allowed gaps std::vector<std::pair<int, int>> staPattern; - for (int gap = 0; gap <= frAlgo.fWData.CurrentIteration()->GetMaxStationGap(); gap++) { + for (int gap = 0; gap <= frWData.CurrentIteration()->GetMaxStationGap(); gap++) { for (int i = 0; i <= gap; i++) { staPattern.push_back(std::make_pair(1 + i, 2 + gap)); } } - const int iStFirst = frAlgo.fWData.CurrentIteration()->GetFirstStationIndex(); + const int iStFirst = frWData.CurrentIteration()->GetFirstStationIndex(); for (int istal = frAlgo.GetParameters().GetNstationsActive() - 2; istal >= iStFirst; istal--) { // start with downstream chambers - const auto& grid = frAlgo.fWData.Grid(istal); + const auto& grid = frWData.Grid(istal); Tindex nGridEntriesL = grid.GetEntries().size(); for (Tindex iel = 0; iel < nGridEntriesL; ++iel) { ca::HitIndex_t ihitl = grid.GetEntries()[iel].GetObjectId(); @@ -358,7 +360,7 @@ void TrackFinderWindow::CaTrackFinderSlice() if (nNeighbours > 0) { assert((int) neighStation >= istal + 1 - && (int) neighStation <= istal + 1 + frAlgo.fWData.CurrentIteration()->GetMaxStationGap()); + && (int) neighStation <= istal + 1 + frWData.CurrentIteration()->GetMaxStationGap()); } unsigned char level = 0; @@ -400,7 +402,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 (frAlgo.fWData.CurrentIteration()->GetTrackFromTripletsFlag()) { + if (frWData.CurrentIteration()->GetTrackFromTripletsFlag()) { min_level = 0; } @@ -425,7 +427,7 @@ void TrackFinderWindow::CaTrackFinderSlice() fvTrackCandidates.clear(); - for (const auto& h : frAlgo.fWData.Hits()) { + for (const auto& h : frWData.Hits()) { frAlgo.fHitKeyToTrack[h.FrontKey()] = -1; frAlgo.fHitKeyToTrack[h.BackKey()] = -1; } @@ -440,19 +442,19 @@ void TrackFinderWindow::CaTrackFinderSlice() for (Tindex itrip = 0; itrip < (Tindex) fvTriplets[istaF].size(); ++itrip) { ca::Triplet& first_trip = (fvTriplets[istaF][itrip]); - const auto& fstTripLHit = frAlgo.fWData.Hit(first_trip.GetLHit()); + const auto& fstTripLHit = frWData.Hit(first_trip.GetLHit()); if (frAlgo.fvHitKeyFlags[fstTripLHit.FrontKey()] || frAlgo.fvHitKeyFlags[fstTripLHit.BackKey()]) { continue; } // skip track candidates that are too short - int minNhits = frAlgo.fWData.CurrentIteration()->GetMinNhits(); + int minNhits = frWData.CurrentIteration()->GetMinNhits(); if (fstTripLHit.Station() == 0) { - minNhits = frAlgo.fWData.CurrentIteration()->GetMinNhitsStation0(); + minNhits = frWData.CurrentIteration()->GetMinNhitsStation0(); } - if (frAlgo.fWData.CurrentIteration()->GetTrackFromTripletsFlag()) { + if (frWData.CurrentIteration()->GetTrackFromTripletsFlag()) { minNhits = 0; } @@ -497,9 +499,9 @@ void TrackFinderWindow::CaTrackFinderSlice() best_tr.SetChi2(best_tr.Chi2() / ndf); if (frAlgo.GetParameters().GetGhostSuppression()) { if (3 == best_tr.NofHits()) { - if (!frAlgo.fWData.CurrentIteration()->GetPrimaryFlag() && (istaF != 0)) + if (!frWData.CurrentIteration()->GetPrimaryFlag() && (istaF != 0)) continue; // too /*short*/ non-MAPS track - if (frAlgo.fWData.CurrentIteration()->GetPrimaryFlag() && (best_tr.Chi2() > 5.0)) continue; + if (frWData.CurrentIteration()->GetPrimaryFlag() && (best_tr.Chi2() > 5.0)) continue; } } fvTrackCandidates.push_back(best_tr); @@ -625,9 +627,9 @@ void TrackFinderWindow::CaTrackFinderSlice() if (!tr.IsAlive()) continue; - if (frAlgo.fWData.CurrentIteration()->GetExtendTracksFlag()) { + if (frWData.CurrentIteration()->GetExtendTracksFlag()) { if (tr.NofHits() < frAlgo.GetParameters().GetNstationsActive()) { - frAlgo.fTrackExtender.ExtendBranch(tr); + fTrackExtender.ExtendBranch(tr); } } @@ -639,11 +641,11 @@ void TrackFinderWindow::CaTrackFinderSlice() frAlgo.fvHitKeyFlags[hit.FrontKey()] = 1; frAlgo.fvHitKeyFlags[hit.BackKey()] = 1; - frAlgo.fWData.RecoHitIndices().push_back(iHit); + frWData.RecoHitIndices().push_back(iHit); } Track t; t.fNofHits = tr.NofHits(); - frAlgo.fWData.RecoTracks().push_back(t); + frWData.RecoTracks().push_back(t); if (0) { // SG debug std::stringstream s; s << "store track " << iCandidate << " chi2= " << tr.Chi2() << "\n"; @@ -660,9 +662,9 @@ void TrackFinderWindow::CaTrackFinderSlice() // suppress strips of suppressed hits - for (unsigned int ih = 0; ih < frAlgo.fWData.Hits().size(); ih++) { - if (frAlgo.fWData.IsHitSuppressed(ih)) { - const ca::Hit& hit = frAlgo.fWData.Hit(ih); + for (unsigned int ih = 0; ih < frWData.Hits().size(); ih++) { + if (frWData.IsHitSuppressed(ih)) { + const ca::Hit& hit = frWData.Hit(ih); frAlgo.fvHitKeyFlags[hit.FrontKey()] = 1; frAlgo.fvHitKeyFlags[hit.BackKey()] = 1; } @@ -670,10 +672,16 @@ void TrackFinderWindow::CaTrackFinderSlice() } // ---- Loop over Track Finder iterations: END -----------------------------------------------------------// // Fit tracks - frAlgo.fTrackFitter.FitCaTracks(); + frAlgo.fMonitorData.StartTimer(ETimer::TrackFitter); + fTrackFitter.FitCaTracks(); + frAlgo.fMonitorData.StopTimer(ETimer::TrackFitter); // Merge clones - frAlgo.fCloneMerger.Exec(frAlgo.fWData.RecoTracks(), frAlgo.fWData.RecoHitIndices()); + fCloneMerger.Exec(frWData.RecoTracks(), frWData.RecoHitIndices()); + + frAlgo.fMonitorData.StartTimer(ETimer::TrackFitter); + fTrackFitter.FitCaTracks(); + frAlgo.fMonitorData.StopTimer(ETimer::TrackFitter); } @@ -697,8 +705,8 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri { // -- finish with current track // add rest of hits - const auto& hitM = frAlgo.fWData.Hit(curr_trip->GetMHit()); - const auto& hitR = frAlgo.fWData.Hit(curr_trip->GetRHit()); + const auto& hitM = frWData.Hit(curr_trip->GetMHit()); + const auto& hitR = frWData.Hit(curr_trip->GetRHit()); if (!(frAlgo.fvHitKeyFlags[hitM.FrontKey()] || frAlgo.fvHitKeyFlags[hitM.BackKey()])) { curr_tr.AddHit(hitM.Id()); @@ -718,7 +726,7 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri } - if (curr_tr.Chi2() > frAlgo.fWData.CurrentIteration()->GetTrackChi2Cut() * ndf) { + if (curr_tr.Chi2() > frWData.CurrentIteration()->GetTrackChi2Cut() * ndf) { return; } @@ -746,7 +754,7 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri fscal dchi2 = 0.; if (!checkTripletMatch(*curr_trip, new_trip, dchi2)) continue; - const auto& hitL = frAlgo.fWData.Hit(new_trip.GetLHit()); // left hit of new triplet + const auto& hitL = frWData.Hit(new_trip.GetLHit()); // left hit of new triplet if (frAlgo.fvHitKeyFlags[hitL.FrontKey()] || frAlgo.fvHitKeyFlags[hitL.BackKey()]) { //hits are used // no used hits allowed -> compare and store track if ((curr_tr.NofHits() > best_tr.NofHits()) @@ -788,7 +796,7 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri ndf = new_L; // TODO: SG: 2 * (new_L + 2) - 5 } - if (new_chi2 > frAlgo.fWData.CurrentIteration()->GetTrackChi2Cut() * ndf) continue; + if (new_chi2 > frWData.CurrentIteration()->GetTrackChi2Cut() * ndf) continue; // add new hit new_tr[ista] = curr_tr; diff --git a/algo/ca/core/tracking/CaTrackFinderWindow.h b/algo/ca/core/tracking/CaTrackFinderWindow.h index 830e3c4685349018d4ae16a23de7d7c41bbbd5df..2813e5ccc22e73344a0a03a3a27453ad823f71ee 100644 --- a/algo/ca/core/tracking/CaTrackFinderWindow.h +++ b/algo/ca/core/tracking/CaTrackFinderWindow.h @@ -10,11 +10,15 @@ #pragma once // include this header only once per compilation unit #include "CaBranch.h" +#include "CaCloneMerger.h" #include "CaGrid.h" #include "CaHit.h" #include "CaSimd.h" +#include "CaTrackExtender.h" +#include "CaTrackFitter.h" #include "CaTrackParam.h" #include "CaVector.h" +#include "CaWindowData.h" namespace cbm::algo::ca { @@ -33,16 +37,16 @@ namespace cbm::algo::ca class TrackFinderWindow { public: /// Default constructor - TrackFinderWindow(ca::Framework& algo); + TrackFinderWindow(ca::Framework& algo, ca::WindowData& wData); /// Destructor ~TrackFinderWindow(); /// Copy constructor - TrackFinderWindow(const TrackFinderWindow&) = delete; + TrackFinderWindow(const TrackFinderWindow&) = default; /// Move constructor - TrackFinderWindow(TrackFinderWindow&&) = delete; + TrackFinderWindow(TrackFinderWindow&&) = default; /// Copy assignment operator TrackFinderWindow& operator=(const TrackFinderWindow&) = delete; @@ -79,6 +83,11 @@ namespace cbm::algo::ca ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class static constexpr bool fDebug = false; // print debug info + + WindowData& frWData; + TrackExtender fTrackExtender; ///< Object of the track extender algorithm + CloneMerger fCloneMerger; ///< Object of the clone merger algorithm + TrackFitter fTrackFitter; ///< Object of the track extender algorithm }; } // namespace cbm::algo::ca diff --git a/algo/ca/core/tracking/CaTrackFitter.cxx b/algo/ca/core/tracking/CaTrackFitter.cxx index d25098a1ecf0619cb2229211ee24bd66b7fa3f1e..26a0f22a1577b97abc1dda35e9ed36d3a623dfc3 100644 --- a/algo/ca/core/tracking/CaTrackFitter.cxx +++ b/algo/ca/core/tracking/CaTrackFitter.cxx @@ -20,7 +20,7 @@ using namespace cbm::algo; // --------------------------------------------------------------------------------------------------------------------- // -TrackFitter::TrackFitter(ca::Framework& algo) : frAlgo(algo) {} +TrackFitter::TrackFitter(ca::Framework& algo, WindowData& wData) : frAlgo(algo), frWData(wData) {} // --------------------------------------------------------------------------------------------------------------------- // @@ -31,7 +31,7 @@ TrackFitter::~TrackFitter() {} void TrackFitter::FitCaTracks() { // LOG(info) << " Start CA Track Fitter "; - int start_hit = 0; // for interation in frAlgo.fWData.RecoHitIndices() + int start_hit = 0; // for interation in frWData.RecoHitIndices() // static ca::FieldValue fldB0, fldB1, fldB2 _fvecalignment; // static ca::FieldRegion fld _fvecalignment; @@ -102,19 +102,19 @@ void TrackFitter::FitCaTracks() mxy[ista].SetCov(1., 0., 1.); } - unsigned short N_vTracks = frAlgo.fWData.RecoTracks().size(); // number of tracks processed per one SSE register + unsigned short N_vTracks = frWData.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] = &frAlgo.fWData.RecoTrack(itrack + i); + t[i] = &frWData.RecoTrack(itrack + i); } // fill the rest of the SIMD vectors with something reasonable for (uint i = nTracks_SIMD; i < fvec::size(); i++) { - t[i] = &frAlgo.fWData.RecoTrack(itrack); + t[i] = &frWData.RecoTrack(itrack); } // get hits of current track @@ -133,7 +133,7 @@ void TrackFitter::FitCaTracks() for (int ih = 0; ih < nHitsTrack; ih++) { - const ca::Hit& hit = frAlgo.GetInputData().GetHit(frAlgo.fWData.RecoHitIndex(start_hit++)); + const ca::Hit& hit = frAlgo.GetInputData().GetHit(frWData.RecoHitIndex(start_hit++)); const int ista = hit.Station(); //if (sta[ista].fieldStatus) { isFieldPresent[iVec] = true; } @@ -279,7 +279,7 @@ void TrackFitter::FitCaTracks() { fitpv.SetMask(fmask::One()); - ca::MeasurementXy<fvec> vtxInfo = frAlgo.fWData.TargetMeasurement(); + ca::MeasurementXy<fvec> vtxInfo = frWData.TargetMeasurement(); vtxInfo.SetDx2(1.e-8); vtxInfo.SetDxy(0.); vtxInfo.SetDy2(1.e-8); diff --git a/algo/ca/core/tracking/CaTrackFitter.h b/algo/ca/core/tracking/CaTrackFitter.h index 75ecd3f186476a48d84046f462fdbbfa622ecab3..2b451763b34267e9d92ad8c480f97daa87a6c6bd 100644 --- a/algo/ca/core/tracking/CaTrackFitter.h +++ b/algo/ca/core/tracking/CaTrackFitter.h @@ -13,6 +13,7 @@ #include "CaSimd.h" #include "CaTrackParam.h" #include "CaVector.h" +#include "CaWindowData.h" namespace cbm::algo::ca @@ -31,16 +32,16 @@ namespace cbm::algo::ca class TrackFitter { public: /// Default constructor - TrackFitter(ca::Framework& algo); + TrackFitter(ca::Framework& algo, WindowData& wData); /// Destructor ~TrackFitter(); /// Copy constructor - TrackFitter(const TrackFitter&) = delete; + TrackFitter(const TrackFitter&) = default; /// Move constructor - TrackFitter(TrackFitter&&) = delete; + TrackFitter(TrackFitter&&) = default; /// Copy assignment operator TrackFitter& operator=(const TrackFitter&) = delete; @@ -61,6 +62,7 @@ namespace cbm::algo::ca /// Data members ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class + ca::WindowData& frWData; }; } // namespace cbm::algo::ca diff --git a/algo/ca/core/tracking/CaTripletConstructor.cxx b/algo/ca/core/tracking/CaTripletConstructor.cxx index 031f0297a1295700538681e08cd2713bca6abc40..abb6514e5fbcf410caa757a2cf62daa70fc7a558 100644 --- a/algo/ca/core/tracking/CaTripletConstructor.cxx +++ b/algo/ca/core/tracking/CaTripletConstructor.cxx @@ -19,7 +19,7 @@ using namespace cbm::algo::ca; -TripletConstructor::TripletConstructor(ca::Framework* algo) { fAlgo = algo; } +TripletConstructor::TripletConstructor(ca::Framework& algo, WindowData& wData) : frAlgo(algo), frWData(wData) {} void TripletConstructor::InitStations(int istal, int istam, int istar) { @@ -29,29 +29,29 @@ void TripletConstructor::InitStations(int istal, int istam, int istar) fIstaM = istam; fIstaR = istar; - if (fIstaM >= fAlgo->fParameters.GetNstationsActive()) { + if (fIstaM >= frAlgo.fParameters.GetNstationsActive()) { return; } - fStaL = &fAlgo->fParameters.GetStation(fIstaL); - fStaM = &fAlgo->fParameters.GetStation(fIstaM); - fStaR = &fAlgo->fParameters.GetStation(fIstaR); + fStaL = &frAlgo.fParameters.GetStation(fIstaL); + fStaM = &frAlgo.fParameters.GetStation(fIstaM); + fStaR = &frAlgo.fParameters.GetStation(fIstaR); { // two stations for approximating the field between the target and the left hit int sta1 = fIstaL; - if (sta1 >= fAlgo->fNfieldStations) { - sta1 = fAlgo->fNfieldStations - 1; + if (sta1 >= frAlgo.fNfieldStations) { + sta1 = frAlgo.fNfieldStations - 1; } if (sta1 < 1) { sta1 = 1; } int sta0 = sta1 / 2; // station between fIstaL and the target - assert(0 <= sta0 && sta0 < sta1 && sta1 < fAlgo->fParameters.GetNstationsActive()); + assert(0 <= sta0 && sta0 < sta1 && sta1 < frAlgo.fParameters.GetNstationsActive()); - fFld0Sta[0] = &fAlgo->fParameters.GetStation(sta0); - fFld0Sta[1] = &fAlgo->fParameters.GetStation(sta1); + fFld0Sta[0] = &frAlgo.fParameters.GetStation(sta0); + fFld0Sta[1] = &frAlgo.fParameters.GetStation(sta1); } { // three stations for approximating the field between the left and the right hit @@ -59,7 +59,7 @@ void TripletConstructor::InitStations(int istal, int istam, int istar) int sta0 = fIstaL; int sta1 = fIstaM; int sta2 = fIstaM + 1; - if (sta2 >= fAlgo->fParameters.GetNstationsActive()) { + if (sta2 >= frAlgo.fParameters.GetNstationsActive()) { sta2 = fIstaM; sta1 = fIstaM - 1; if (sta0 >= sta1) { @@ -69,17 +69,17 @@ void TripletConstructor::InitStations(int istal, int istam, int istar) sta0 = 0; } } - if (fAlgo->fParameters.GetNstationsActive() >= 3) { - assert(0 <= sta0 && sta0 < sta1 && sta1 < sta2 && sta2 < fAlgo->fParameters.GetNstationsActive()); + if (frAlgo.fParameters.GetNstationsActive() >= 3) { + assert(0 <= sta0 && sta0 < sta1 && sta1 < sta2 && sta2 < frAlgo.fParameters.GetNstationsActive()); } - fFld1Sta[0] = &fAlgo->fParameters.GetStation(sta0); - fFld1Sta[1] = &fAlgo->fParameters.GetStation(sta1); - fFld1Sta[2] = &fAlgo->fParameters.GetStation(sta2); + fFld1Sta[0] = &frAlgo.fParameters.GetStation(sta0); + fFld1Sta[1] = &frAlgo.fParameters.GetStation(sta1); + fFld1Sta[2] = &frAlgo.fParameters.GetStation(sta2); } - fFit.SetParticleMass(fAlgo->GetDefaultParticleMass()); - if (fAlgo->fWData.CurrentIteration()->GetElectronFlag()) { + fFit.SetParticleMass(frAlgo.GetDefaultParticleMass()); + if (frWData.CurrentIteration()->GetElectronFlag()) { fFit.SetParticleMass(constants::phys::ElectronMass); } fFit.SetMask(fmask::One()); @@ -96,7 +96,7 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c if (!fIsInitialized) return; fIhitL = iHitL; - const auto& hitl = fAlgo->fWData.Hit(fIhitL); + const auto& hitl = frWData.Hit(fIhitL); // fit a straight line through the target and the left hit. { @@ -104,7 +104,7 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c /// Get the field approximation. Add the target to parameters estimation. /// Propagaete to the middle station of a triplet. - fFit.SetQp0(fAlgo->fWData.CurrentIteration()->GetMaxQp()); + fFit.SetQp0(frWData.CurrentIteration()->GetMaxQp()); ca::FieldValue<fvec> lB, mB, cB, rB _fvecalignment; ca::FieldValue<fvec> l_B_global, centB_global _fvecalignment; @@ -117,10 +117,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 - fAlgo->fWData.TargZ()); + const fvec dzli = 1.f / (zl - frAlgo.GetParameters().GetTargetPositionZ()); - const fvec tx = (xl - fAlgo->fWData.TargX()) * dzli; - const fvec ty = (yl - fAlgo->fWData.TargY()) * dzli; + const fvec tx = (xl - frAlgo.GetParameters().GetTargetPositionX()) * dzli; + const fvec ty = (yl - frAlgo.GetParameters().GetTargetPositionY()) * dzli; TrackParamV& T = fFit.Tr(); @@ -133,20 +133,20 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c T.Time() = time; T.Vi() = constants::phys::SpeedOfLightInv; - const fvec maxSlopePV = fAlgo->fWData.CurrentIteration()->GetMaxSlopePV(); - const fvec maxQp = fAlgo->fWData.CurrentIteration()->GetMaxQp(); + const fvec maxSlopePV = frWData.CurrentIteration()->GetMaxSlopePV(); + const fvec maxQp = frWData.CurrentIteration()->GetMaxQp(); fvec txErr2 = maxSlopePV * maxSlopePV / fvec(9.); fvec qpErr2 = maxQp * maxQp / fvec(9.); T.ResetErrors(1., 1., txErr2, txErr2, qpErr2, (staL().timeInfo ? timeEr2 : 1.e6), 1.e10); - T.InitVelocityRange(1. / fAlgo->fWData.CurrentIteration()->GetMaxQp()); + T.InitVelocityRange(1. / frWData.CurrentIteration()->GetMaxQp()); T.C00() = hitl.dX2(); T.C10() = hitl.dXY(); T.C11() = hitl.dY2(); - T.Ndf() = (fAlgo->fWData.CurrentIteration()->GetPrimaryFlag()) ? fvec(2.) : fvec(0.); + T.Ndf() = (frWData.CurrentIteration()->GetPrimaryFlag()) ? fvec(2.) : fvec(0.); T.NdfTime() = (staL().timeInfo ? 0 : -1); // NDF = number of track parameters (6: x, y, tx, ty, qp, time) - number of measured parameters (3: x, y, time) on station or (2: x, y) on target @@ -157,7 +157,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(fAlgo->fWData.TargB(), fAlgo->fWData.TargZ(), B0, fFld0Sta[0]->fZ, B1, fFld0Sta[1]->fZ); + fld0.Set(frWData.TargB(), frAlgo.GetParameters().GetTargetPositionZ(), B0, fFld0Sta[0]->fZ, B1, fFld0Sta[1]->fZ); } { // field, made by the left hit, the middle station and the right station @@ -170,9 +170,9 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c // add the target constraint - fFit.FilterWithTargetAtLine(fAlgo->fWData.TargZ(), fAlgo->fWData.TargetMeasurement(), fld0); + fFit.FilterWithTargetAtLine(frAlgo.GetParameters().GetTargetPositionZ(), frWData.TargetMeasurement(), fld0); - fFit.MultipleScattering(fAlgo->fParameters.GetMaterialThickness(fIstaL, T.GetX(), T.GetY())); + fFit.MultipleScattering(frAlgo.fParameters.GetMaterialThickness(fIstaL, T.GetX(), T.GetY())); // extrapolate to the middle hit @@ -192,15 +192,15 @@ void TripletConstructor::FindDoublets() /// Find the doublets. Reformat data into portions of doublets. int iMC = -1; - if (fAlgo->fParameters.DevIsMatchDoubletsViaMc()) { - iMC = fAlgo->GetMcTrackIdForWindowHit(fIhitL); + if (frAlgo.fParameters.DevIsMatchDoubletsViaMc()) { + iMC = frAlgo.GetMcTrackIdForWindowHit(fIhitL); if (iMC < 0) { return; } } - CollectHits(fTrL, fIstaM, fAlgo->fWData.CurrentIteration()->GetDoubletChi2Cut(), iMC, fHitsM_2, - fAlgo->fParameters.GetMaxDoubletsPerSinglet()); + CollectHits(fTrL, fIstaM, frWData.CurrentIteration()->GetDoubletChi2Cut(), iMC, fHitsM_2, + frAlgo.fParameters.GetMaxDoubletsPerSinglet()); } @@ -215,15 +215,15 @@ void TripletConstructor::FitDoublets() fTracks_2.clear(); fTracks_2.reserve(hitsMtmp.size()); - bool isMomentumFitted = (fAlgo->fIsTargetField || (staL().fieldStatus != 0) || (staM().fieldStatus != 0)); + bool isMomentumFitted = (frAlgo.fIsTargetField || (staL().fieldStatus != 0) || (staM().fieldStatus != 0)); - double maxQp = fAlgo->fWData.CurrentIteration()->GetMaxQp(); + double maxQp = frWData.CurrentIteration()->GetMaxQp(); for (unsigned int i2 = 0; i2 < hitsMtmp.size(); i2++) { const ca::HitIndex_t indM = hitsMtmp[i2]; - const ca::Hit& hitm = fAlgo->fWData.Hit(indM); + const ca::Hit& hitm = frWData.Hit(indM); - if (fAlgo->fWData.IsHitSuppressed(indM)) { + if (frWData.IsHitSuppressed(indM)) { continue; } @@ -246,12 +246,12 @@ void TripletConstructor::FitDoublets() fFit.SetQp0(isMomentumFitted ? fFit.Tr().GetQp() : maxQp); - fFit.MultipleScattering(fAlgo->fParameters.GetMaterialThickness(fIstaM, T2.GetX(), T2.Y())); + fFit.MultipleScattering(frAlgo.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 != fAlgo->fTrackingMode) { + if (ca::Framework::kMcbm != frAlgo.fTrackingMode) { // TODO: SG: adjust cuts, put them to parameters fscal tx = T2.Tx()[0]; @@ -261,7 +261,7 @@ void TripletConstructor::FitDoublets() for (unsigned int iClone = i2 + 1; iClone < hitsMtmp.size(); iClone++) { int indClone = hitsMtmp[iClone]; - const ca::Hit& hitClone = fAlgo->fWData.Hit(indClone); + const ca::Hit& hitClone = frWData.Hit(indClone); fscal dz = hitClone.Z() - T2.Z()[0]; @@ -282,13 +282,13 @@ void TripletConstructor::FitDoublets() continue; } - if (fAlgo->fParameters.DevIsSuppressOverlapHitsViaMc()) { - int iMC = fAlgo->GetMcTrackIdForWindowHit(fIhitL); - if ((iMC != fAlgo->GetMcTrackIdForWindowHit(indM)) || (iMC != fAlgo->GetMcTrackIdForWindowHit(indClone))) { + if (frAlgo.fParameters.DevIsSuppressOverlapHitsViaMc()) { + int iMC = frAlgo.GetMcTrackIdForWindowHit(fIhitL); + if ((iMC != frAlgo.GetMcTrackIdForWindowHit(indM)) || (iMC != frAlgo.GetMcTrackIdForWindowHit(indClone))) { continue; } } - fAlgo->fWData.SuppressHit(indClone); + frWData.SuppressHit(indClone); } } @@ -300,7 +300,7 @@ void TripletConstructor::FitDoublets() void TripletConstructor::FindTriplets() { - if (fIstaR >= fAlgo->fParameters.GetNstationsActive()) { + if (fIstaR >= frAlgo.fParameters.GetNstationsActive()) { return; } FindRightHit(); @@ -317,10 +317,10 @@ void TripletConstructor::FindRightHit() fFit.SetQp0(fvec(0.)); - if (fIstaR >= fAlgo->fParameters.GetNstationsActive()) return; + if (fIstaR >= frAlgo.fParameters.GetNstationsActive()) return; { - int maxTriplets = fHitsM_2.size() * fAlgo->fParameters.GetMaxTripletPerDoublets(); + int maxTriplets = fHitsM_2.size() * frAlgo.fParameters.GetMaxTripletPerDoublets(); fHitsM_3.clear(); fHitsR_3.clear(); fHitsM_3.reserve(maxTriplets); @@ -328,8 +328,8 @@ void TripletConstructor::FindRightHit() } // ---- Add the middle hits to parameters estimation. Propagate to right station. ---- - double maxSlope = fAlgo->fWData.CurrentIteration()->GetMaxSlope(); - double tripletChi2Cut = fAlgo->fWData.CurrentIteration()->GetTripletChi2Cut(); + double maxSlope = frWData.CurrentIteration()->GetMaxSlope(); + double tripletChi2Cut = frWData.CurrentIteration()->GetTripletChi2Cut(); for (unsigned int i2 = 0; i2 < fHitsM_2.size(); i2++) { fFit.SetTrack(fTracks_2[i2]); @@ -341,15 +341,15 @@ void TripletConstructor::FindRightHit() if constexpr (fDebugDublets) { ca::HitIndex_t iwhit[2] = {fIhitL, fHitsM_2[i2]}; - ca::HitIndex_t ihit[2] = {fAlgo->fWData.Hit(iwhit[0]).Id(), fAlgo->fWData.Hit(iwhit[1]).Id()}; + ca::HitIndex_t ihit[2] = {frWData.Hit(iwhit[0]).Id(), frWData.Hit(iwhit[1]).Id()}; int ih0 = ihit[0]; int ih1 = ihit[1]; - const ca::Hit& h0 = fAlgo->fInputData.GetHit(ih0); - const ca::Hit& h1 = fAlgo->fInputData.GetHit(ih1); + const ca::Hit& h0 = frAlgo.fInputData.GetHit(ih0); + const ca::Hit& h1 = frAlgo.fInputData.GetHit(ih1); LOG(info) << "\n====== Extrapolated Doublet : " - << " iter " << fAlgo->fCurrentIterationIndex << " hits: {" << fIstaL << "/" << ih0 << " " << fIstaM + << " iter " << frAlgo.fCurrentIterationIndex << " hits: {" << fIstaL << "/" << ih0 << " " << fIstaM << "/" << ih1 << " " << fIstaR << "/?} xyz: {" << h0.X() << " " << h0.Y() << " " << h0.Z() << "}, {" << h1.X() << " " << h1.Y() << " " << h1.Z() << "} chi2 " << T2.GetChiSq()[0] << " ndf " << T2.Ndf()[0] << " chi2time " << T2.ChiSqTime()[0] << " ndfTime " << T2.NdfTime()[0]; @@ -360,7 +360,7 @@ void TripletConstructor::FindRightHit() bool isDoubletGood = false; int iMC = -1; do { - if (fAlgo->kSts == fAlgo->fTrackingMode && (T2.C44()[0] < 0)) { + if (frAlgo.kSts == frAlgo.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 +372,10 @@ void TripletConstructor::FindRightHit() if (fabs(T2.Ty()[0]) > maxSlope) { break; } - if (fAlgo->fParameters.DevIsMatchTripletsViaMc()) { + if (frAlgo.fParameters.DevIsMatchTripletsViaMc()) { int indM = fHitsM_2[i2]; - iMC = fAlgo->GetMcTrackIdForWindowHit(fIhitL); - if (iMC < 0 || iMC != fAlgo->GetMcTrackIdForWindowHit(indM)) { + iMC = frAlgo.GetMcTrackIdForWindowHit(fIhitL); + if (iMC < 0 || iMC != frAlgo.GetMcTrackIdForWindowHit(indM)) { break; } } @@ -397,10 +397,10 @@ void TripletConstructor::FindRightHit() } Vector<ca::HitIndex_t> collectedHits; - CollectHits(T2, fIstaR, tripletChi2Cut, iMC, collectedHits, fAlgo->fParameters.GetMaxTripletPerDoublets()); + CollectHits(T2, fIstaR, tripletChi2Cut, iMC, collectedHits, frAlgo.fParameters.GetMaxTripletPerDoublets()); - if (collectedHits.size() >= fAlgo->fParameters.GetMaxTripletPerDoublets()) { - LOG(debug) << "Ca: GetMaxTripletPerDoublets==" << fAlgo->fParameters.GetMaxTripletPerDoublets() + if (collectedHits.size() >= frAlgo.fParameters.GetMaxTripletPerDoublets()) { + LOG(debug) << "Ca: GetMaxTripletPerDoublets==" << frAlgo.fParameters.GetMaxTripletPerDoublets() << " reached in findTripletsStep0()"; // reject already created triplets for this doublet collectedHits.clear(); @@ -411,11 +411,11 @@ void TripletConstructor::FindRightHit() for (unsigned int ih = 0; ih < collectedHits.size(); ih++) { ca::HitIndex_t irh = collectedHits[ih]; if constexpr (fDebugDublets) { - ca::HitIndex_t ihit = fAlgo->fWData.Hit(irh).Id(); - const ca::Hit& h = fAlgo->fInputData.GetHit(ihit); + ca::HitIndex_t ihit = frWData.Hit(irh).Id(); + const ca::Hit& h = frAlgo.fInputData.GetHit(ihit); LOG(info) << " hit " << ihit << " " << h.ToString(); } - if (fAlgo->fWData.IsHitSuppressed(irh)) { + if (frWData.IsHitSuppressed(irh)) { if constexpr (fDebugDublets) { LOG(info) << " the hit is suppressed"; } @@ -451,11 +451,11 @@ void TripletConstructor::FitTriplets() ca::TrackFit fit; fit.SetMask(fmask::One()); - if (fAlgo->fWData.CurrentIteration()->GetElectronFlag()) { + if (frWData.CurrentIteration()->GetElectronFlag()) { fit.SetParticleMass(constants::phys::ElectronMass); } else { - fit.SetParticleMass(fAlgo->GetDefaultParticleMass()); + fit.SetParticleMass(frAlgo.GetDefaultParticleMass()); } const int NHits = 3; // triplets @@ -464,7 +464,7 @@ void TripletConstructor::FitTriplets() ca::Station<fvec> sta[3]; for (int is = 0; is < NHits; ++is) { - sta[is] = fAlgo->fParameters.GetStation(ista[is]); + sta[is] = frAlgo.fParameters.GetStation(ista[is]); }; bool isMomentumFitted = ((staL().fieldStatus != 0) || (staM().fieldStatus != 0) || (staR().fieldStatus != 0)); @@ -472,19 +472,18 @@ void TripletConstructor::FitTriplets() fvec ndfTrackModel = 4; // straight line ndfTrackModel += isMomentumFitted ? 1 : 0; // track with momentum - const fvec maxQp = fAlgo->fWData.CurrentIteration()->GetMaxQp(); + const fvec maxQp = frWData.CurrentIteration()->GetMaxQp(); for (int i3 = 0; i3 < n3; ++i3) { // prepare data ca::HitIndex_t iwhit[NHits] = {fIhitL, fHitsM_3[i3], fHitsR_3[i3]}; - ca::HitIndex_t ihit[NHits] = {fAlgo->fWData.Hit(iwhit[0]).Id(), fAlgo->fWData.Hit(iwhit[1]).Id(), - fAlgo->fWData.Hit(iwhit[2]).Id()}; + ca::HitIndex_t ihit[NHits] = {frWData.Hit(iwhit[0]).Id(), frWData.Hit(iwhit[1]).Id(), frWData.Hit(iwhit[2]).Id()}; - if (fAlgo->fParameters.DevIsMatchTripletsViaMc()) { - int mc1 = fAlgo->GetMcTrackIdForCaHit(ihit[0]); - int mc2 = fAlgo->GetMcTrackIdForCaHit(ihit[1]); - int mc3 = fAlgo->GetMcTrackIdForCaHit(ihit[2]); + if (frAlgo.fParameters.DevIsMatchTripletsViaMc()) { + int mc1 = frAlgo.GetMcTrackIdForCaHit(ihit[0]); + int mc2 = frAlgo.GetMcTrackIdForCaHit(ihit[1]); + int mc3 = frAlgo.GetMcTrackIdForCaHit(ihit[2]); if ((mc1 != mc2) || (mc1 != mc3)) { continue; } @@ -495,7 +494,7 @@ void TripletConstructor::FitTriplets() ca::MeasurementXy<fvec> mxy[NHits]; for (int ih = 0; ih < NHits; ++ih) { - const ca::Hit& hit = fAlgo->fInputData.GetHit(ihit[ih]); + const ca::Hit& hit = frAlgo.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(); @@ -519,7 +518,7 @@ void TripletConstructor::FitTriplets() }; fld.Set(B[0], sta[0].fZ, B[1], sta[1].fZ, B[2], sta[2].fZ); - fldTarget.Set(fAlgo->fWData.TargB(), fAlgo->fWData.TargZ(), B[0], sta[0].fZ, B[1], sta[1].fZ); + fldTarget.Set(frWData.TargB(), frAlgo.GetParameters().GetTargetPositionZ(), B[0], sta[0].fZ, B[1], sta[1].fZ); TrackParamV& T = fit.Tr(); @@ -557,11 +556,11 @@ void TripletConstructor::FitTriplets() T.NdfTime() = sta[ih0].timeInfo ? 0 : -1; // add the target constraint - fit.FilterWithTargetAtLine(fAlgo->fWData.TargZ(), fAlgo->fWData.TargetMeasurement(), fldTarget); + fit.FilterWithTargetAtLine(frAlgo.GetParameters().GetTargetPositionZ(), frWData.TargetMeasurement(), fldTarget); for (int ih = 1; ih < NHits; ++ih) { fit.Extrapolate(z[ih], fld); - auto radThick = fAlgo->fParameters.GetMaterialThickness(ista[ih], T.X(), T.Y()); + auto radThick = frAlgo.fParameters.GetMaterialThickness(ista[ih], T.X(), T.Y()); fit.MultipleScattering(radThick); fit.EnergyLossCorrection(radThick, fvec(-1.f)); fit.FilterXY(mxy[ih]); @@ -595,7 +594,7 @@ void TripletConstructor::FitTriplets() for (int ih = NHits - 2; ih >= 0; --ih) { fit.Extrapolate(z[ih], fld); - auto radThick = fAlgo->fParameters.GetMaterialThickness(ista[ih], T.X(), T.Y()); + auto radThick = frAlgo.fParameters.GetMaterialThickness(ista[ih], T.X(), T.Y()); fit.MultipleScattering(radThick); fit.EnergyLossCorrection(radThick, fvec(1.f)); fit.FilterXY(mxy[ih]); @@ -610,23 +609,23 @@ void TripletConstructor::FitTriplets() int ih0 = ihit[0]; int ih1 = ihit[1]; int ih2 = ihit[2]; - int mc1 = fAlgo->GetMcTrackIdForCaHit(ih0); - int mc2 = fAlgo->GetMcTrackIdForCaHit(ih1); - int mc3 = fAlgo->GetMcTrackIdForCaHit(ih2); + int mc1 = frAlgo.GetMcTrackIdForCaHit(ih0); + int mc2 = frAlgo.GetMcTrackIdForCaHit(ih1); + int mc3 = frAlgo.GetMcTrackIdForCaHit(ih2); if (1 || (mc1 >= 0) && (mc1 == mc2) && (mc1 == mc3)) { - const ca::Hit& h0 = fAlgo->fInputData.GetHit(ih0); - const ca::Hit& h1 = fAlgo->fInputData.GetHit(ih1); - const ca::Hit& h2 = fAlgo->fInputData.GetHit(ih2); + const ca::Hit& h0 = frAlgo.fInputData.GetHit(ih0); + const ca::Hit& h1 = frAlgo.fInputData.GetHit(ih1); + const ca::Hit& h2 = frAlgo.fInputData.GetHit(ih2); //const CbmL1MCTrack& mctr = CbmL1::Instance()->GetMcTracks()[mc1]; LOG(info) << "== fitted triplet: " - << " iter " << fAlgo->fCurrentIterationIndex << " hits: {" << fIstaL << "/" << ih0 << " " << fIstaM + << " iter " << frAlgo.fCurrentIterationIndex << " hits: {" << fIstaL << "/" << ih0 << " " << fIstaM << "/" << ih1 << " " << fIstaR << "/" << ih2 << "} xyz: {" << h0.X() << " " << h0.Y() << " " << h0.Z() << "}, {" << h1.X() << " " << h1.Y() << " " << h1.Z() << "}, {" << h2.X() << ", " << h2.Y() << ", " << h2.Z() << "} chi2 " << T.GetChiSq()[0] << " ndf " << T.Ndf()[0] << " chi2time " << T.ChiSqTime()[0] << " ndfTime " << T.NdfTime()[0]; /* cbm::ca::tools::Debugger::Instance().FillNtuple( - "triplets", mctr.iEvent, fAlgo->fCurrentIterationIndex, ih0, h0.X(), h0.Y(), h0.Z(), ih1, h1.X(), h1.Y(), + "triplets", mctr.iEvent, frAlgo.fCurrentIterationIndex, ih0, h0.X(), h0.Y(), h0.Z(), ih1, h1.X(), h1.Y(), h1.Z(), ih2, h2.X(), h2.Y(), h2.Z(), mc1, fIstaL, mctr.p, mctr.x, mctr.y, mctr.z, (float) T.GetChiSq()[0], (float) T.Ndf()[0], (float) T.ChiSqTime()[0], (float) T.NdfTime()[0]); */ @@ -649,7 +648,7 @@ void TripletConstructor::StoreTriplets() fTriplets.clear(); fTriplets.reserve(n3); - const float tripletFinalChi2Cut = fAlgo->fWData.CurrentIteration()->GetTripletFinalChi2Cut(); + const float tripletFinalChi2Cut = frWData.CurrentIteration()->GetTripletFinalChi2Cut(); for (Tindex i3 = 0; i3 < n3; ++i3) { TrackParamV& T3 = fTracks_3[i3]; @@ -662,11 +661,11 @@ void TripletConstructor::StoreTriplets() const ca::HitIndex_t ihitm = fHitsM_3[i3]; const ca::HitIndex_t ihitr = fHitsR_3[i3]; - CBMCA_DEBUG_ASSERT(ihitl < fAlgo->fWData.HitStartIndexOnStation(fIstaL + 1)); - CBMCA_DEBUG_ASSERT(ihitm < fAlgo->fWData.HitStartIndexOnStation(fIstaM + 1)); - CBMCA_DEBUG_ASSERT(ihitr < fAlgo->fWData.HitStartIndexOnStation(fIstaR + 1)); + CBMCA_DEBUG_ASSERT(ihitl < frWData.HitStartIndexOnStation(fIstaL + 1)); + CBMCA_DEBUG_ASSERT(ihitm < frWData.HitStartIndexOnStation(fIstaM + 1)); + CBMCA_DEBUG_ASSERT(ihitr < frWData.HitStartIndexOnStation(fIstaR + 1)); - if (!fAlgo->fWData.CurrentIteration()->GetTrackFromTripletsFlag()) { + if (!frWData.CurrentIteration()->GetTrackFromTripletsFlag()) { if (chi2 > tripletFinalChi2Cut) { continue; } @@ -701,7 +700,7 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons collectedHits.clear(); collectedHits.reserve(maxNhits); - const ca::Station<fvec>& sta = fAlgo->fParameters.GetStation(iSta); + const ca::Station<fvec>& sta = frAlgo.fParameters.GetStation(iSta); fFit.SetTrack(Tr); TrackParamV& T = fFit.Tr(); @@ -713,8 +712,8 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons const fscal timeError2 = T.C55()[0]; const fscal time = T.Time()[0]; - const auto& grid = fAlgo->fWData.Grid(iSta); - const fvec maxDZ = fAlgo->fWData.CurrentIteration()->GetMaxDZ(); + const auto& grid = frWData.Grid(iSta); + const fvec maxDZ = frWData.CurrentIteration()->GetMaxDZ(); ca::GridArea area(grid, T.X()[0], T.Y()[0], (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 constexpr (fDebugCollectHits) { @@ -722,7 +721,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 (fAlgo->fParameters.DevIsIgnoreHitSearchAreas()) { + if (frAlgo.fParameters.DevIsIgnoreHitSearchAreas()) { area.DoLoopOverEntireGrid(); } @@ -731,17 +730,17 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons while (area.GetNextObjectId(ih) && ((int) collectedHits.size() < maxNhits)) { // loop over station hits - if (fAlgo->fWData.IsHitSuppressed(ih)) { + if (frWData.IsHitSuppressed(ih)) { continue; } - const ca::Hit& hit = fAlgo->fWData.Hit(ih); + const ca::Hit& hit = frWData.Hit(ih); if constexpr (fDebugCollectHits) { LOG(info) << "hit in the search area: " << hit.ToString(); LOG(info) << " check the hit.. "; } if (iMC >= 0) { // match via MC - if (iMC != fAlgo->GetMcTrackIdForWindowHit(ih)) { + if (iMC != frAlgo.GetMcTrackIdForWindowHit(ih)) { if constexpr (fDebugCollectHits) { LOG(info) << " hit mc does not match"; } @@ -805,7 +804,7 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons auto [chi2x, chi2u] = ca::TrackFit::GetChi2XChi2U(mxy, x, y, C00, C10, C11); // TODO: adjust the cut, cut on chi2x & chi2u separately - if (!fAlgo->fWData.CurrentIteration()->GetTrackFromTripletsFlag()) { + if (!frWData.CurrentIteration()->GetTrackFromTripletsFlag()) { if (chi2x[0] > chi2Cut) { if constexpr (fDebugCollectHits) { LOG(info) << " hit chi2X is too large"; diff --git a/algo/ca/core/tracking/CaTripletConstructor.h b/algo/ca/core/tracking/CaTripletConstructor.h index 524aab6f42228aa8d2e969f1e80a8bef9b3c1af3..de18a4d718613d726bba6a5148b64b9d71f50e10 100644 --- a/algo/ca/core/tracking/CaTripletConstructor.h +++ b/algo/ca/core/tracking/CaTripletConstructor.h @@ -16,6 +16,7 @@ #include "CaTrackParam.h" #include "CaTriplet.h" #include "CaVector.h" +#include "CaWindowData.h" namespace cbm::algo::ca { @@ -33,7 +34,7 @@ namespace cbm::algo::ca /// Constructor /// \param nThreads Number of threads for multi-threaded mode - TripletConstructor(ca::Framework* algo); + TripletConstructor(ca::Framework& algo, ca::WindowData& wData); /// Copy constructor TripletConstructor(const TripletConstructor&) = delete; @@ -88,8 +89,8 @@ namespace cbm::algo::ca const ca::Station<fvec>& staR() { return *fStaR; } private: - bool fIsInitialized{false}; ///< is initialized; - ca::Framework* fAlgo{nullptr}; ///< pointer to ca::Framework object + 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 @@ -118,6 +119,7 @@ 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/L1Algo/utils/L1AlgoDraw.cxx b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx index 5cf56fb5f334952e173cd2df5a1706a07bc224d2..a203bc8bfcdef5038e04422355c3044e358514dd 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->fWData.Hits().size(); + auto nHits = algo->fvWData[0].Hits().size(); vHits.reserve(nHits); vHitsQa.reserve(nHits); - for (const auto& hit : algo->fWData.Hits()) { + for (const auto& hit : algo->fvWData[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->fWData.HitStartIndexOnStation(i); - HitsStopIndex[i] = HitsStartIndex[i] + algo->fWData.NofHitsOnStation(i); + HitsStartIndex[i] = algo->fvWData[0].HitStartIndexOnStation(i); + HitsStopIndex[i] = HitsStartIndex[i] + algo->fvWData[0].NofHitsOnStation(i); vStations[i] = algo->GetParameters().GetStation(i); } } @@ -245,8 +245,9 @@ void L1AlgoDraw::DrawRecoTracks() // CbmL1 &L1 = *CbmL1::Instance(); int curRecoHit = 0; - cbm::algo::ca::Vector<ca::HitIndex_t>& recoHits = algo->fWData.RecoHitIndices(); - for (vector<Track>::iterator it = algo->fWData.RecoTracks().begin(); it != algo->fWData.RecoTracks().end(); ++it) { + 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; // if (nHits > 5) continue; // draw clones @@ -484,8 +485,8 @@ void L1AlgoDraw::DrawDoublet(int il, int ir) void L1AlgoDraw::DrawInfo() { cout << " vHits.size = " << algo->GetInputData().GetNhits() << endl; - cout << " vRecoHits.size = " << algo->fWData.RecoHitIndices().size() << endl; - cout << " vTracks.size = " << algo->fWData.RecoTracks().size() << endl; + cout << " vRecoHits.size = " << algo->fvWData[0].RecoHitIndices().size() << endl; + cout << " vTracks.size = " << algo->fvWData[0].RecoTracks().size() << endl; } void L1AlgoDraw::DrawTarget() diff --git a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h index 7023d4dc0e69de43984b09a28dbf09de0a621b57..fa0857c561462efa6ad3ae9e2e7aec9eef2cf811 100644 --- a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h +++ b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h @@ -202,13 +202,9 @@ void L1AlgoEfficiencyPerformance<NHits>::FillMC() trlet.mcTrackId = mcTrk.GetId(); trlet.mcMotherTrackId = mcTrk.GetMotherId(); trlet.p = mcTrk.GetP(); - if (fL1->fpAlgo->fWData.CurrentIteration()) { - if (mcTrk.GetP() > 1. / fL1->fpAlgo->fWData.CurrentIteration()->GetMaxQp()) { - trlet.SetAsReconstructable(); - } - } - else { - LOG(error) << "L1AlgoEfficiencyPerformance::FillMC(): attempt to call fpCurrentIteration, but it is a nullptr"; + const fscal maxQp = 1.0 / 0.5; // NOTE: Default value from ca::Iteration + if (mcTrk.GetP() > 1. / maxQp) { + trlet.SetAsReconstructable(); } mcTracklets.push_back(trlet);