From 649a5266612007f8f30d58684be2353ce1a7a8d3 Mon Sep 17 00:00:00 2001 From: Matevz Tadel Date: Tue, 26 Oct 2021 10:03:11 -0700 Subject: [PATCH] For clone-engine keep candidate state at last valid hit. Add option to export candidates without missing hit information, use it for CMSSW runOneIteration. In backward fit avoid re-application of material effects when a track does not have a valid hit on current layer. --- mkFit/HitStructures.cc | 20 +++++++++++------- mkFit/HitStructures.h | 2 +- mkFit/MkBase.h | 10 +++++---- mkFit/MkBuilder.cc | 13 +++++++----- mkFit/MkBuilder.h | 4 ++-- mkFit/MkFinder.cc | 21 +++++++++++++------ mkFit/PropagationMPlex.cc | 44 +++++++++++++++++++++++---------------- mkFit/PropagationMPlex.h | 6 ++++-- mkFit/buildtestMPlex.cc | 2 +- 9 files changed, 76 insertions(+), 46 deletions(-) diff --git a/mkFit/HitStructures.cc b/mkFit/HitStructures.cc index beb37612..03b95375 100644 --- a/mkFit/HitStructures.cc +++ b/mkFit/HitStructures.cc @@ -445,27 +445,33 @@ EventOfHits::EventOfHits(const TrackerInfo &trk_inf) : // TrackCand //============================================================================== -Track TrackCand::exportTrack() const +Track TrackCand::exportTrack(bool remove_missing_hits) const { // printf("TrackCand::exportTrack label=%5d, total_hits=%2d, overlaps=%2d -- n_seed_hits=%d,prod_type=%d\n", // label(), nTotalHits(), nOverlapHits_, getNSeedHits(), (int)prodType()); Track res(*this); - res.resizeHits(nTotalHits(), nFoundHits()); + res.resizeHits(remove_missing_hits ? nFoundHits() : nTotalHits(), nFoundHits()); res.setNOverlapHits(nOverlapHits()); int nh = nTotalHits(); int ch = lastHitIdx_; + int good_hits_pos = nFoundHits(); while (--nh >= 0) { HoTNode& hot_node = m_comb_candidate->m_hots[ch]; - - res.setHitIdxAtPos(nh, hot_node.m_hot); - + if (remove_missing_hits) + { + if (hot_node.m_hot.index >= 0) + res.setHitIdxAtPos(--good_hits_pos, hot_node.m_hot); + } + else + { + res.setHitIdxAtPos(nh, hot_node.m_hot); + } // printf(" nh=%2d, ch=%d, idx=%d lyr=%d prev_idx=%d\n", // nh, ch, hot_node.m_hot.index, hot_node.m_hot.layer, hot_node.m_prev_idx); - - ch = hot_node.m_prev_idx; + ch = hot_node.m_prev_idx; } return res; diff --git a/mkFit/HitStructures.h b/mkFit/HitStructures.h index 6e8e2c95..06ba640a 100644 --- a/mkFit/HitStructures.h +++ b/mkFit/HitStructures.h @@ -528,7 +528,7 @@ class TrackCand : public TrackBase void incOverlapCount() { ++nOverlapHits_; } - Track exportTrack() const; + Track exportTrack(bool remove_missing_hits=false) const; void resetShortTrack() { score_ = getScoreWorstPossible(); m_comb_candidate = nullptr; } diff --git a/mkFit/MkBase.h b/mkFit/MkBase.h index 7bebd4a5..59560619 100644 --- a/mkFit/MkBase.h +++ b/mkFit/MkBase.h @@ -45,7 +45,8 @@ class MkBase Err[iP], Par[iP], N_proc, pf); } - void PropagateTracksToHitR(const MPlexHV& par, const int N_proc, const PropagationFlags pf) + void PropagateTracksToHitR(const MPlexHV& par, const int N_proc, const PropagationFlags pf, + const MPlexQI *noMatEffPtr=nullptr) { MPlexQF msRad; #pragma omp simd @@ -55,7 +56,7 @@ class MkBase } propagateHelixToRMPlex(Err[iC], Par[iC], Chg, msRad, - Err[iP], Par[iP], N_proc, pf); + Err[iP], Par[iP], N_proc, pf, noMatEffPtr); } //---------------------------------------------------------------------------- @@ -73,7 +74,8 @@ class MkBase Err[iP], Par[iP], N_proc, pf); } - void PropagateTracksToHitZ(const MPlexHV& par, const int N_proc, const PropagationFlags pf) + void PropagateTracksToHitZ(const MPlexHV& par, const int N_proc, const PropagationFlags pf, + const MPlexQI *noMatEffPtr=nullptr) { MPlexQF msZ; #pragma omp simd @@ -83,7 +85,7 @@ class MkBase } propagateHelixToZMPlex(Err[iC], Par[iC], Chg, msZ, - Err[iP], Par[iP], N_proc, pf); + Err[iP], Par[iP], N_proc, pf, noMatEffPtr); } void PropagateTracksToPCAZ(const int N_proc, const PropagationFlags pf) diff --git a/mkFit/MkBuilder.cc b/mkFit/MkBuilder.cc index eef2cfbf..d2f437f6 100644 --- a/mkFit/MkBuilder.cc +++ b/mkFit/MkBuilder.cc @@ -303,14 +303,14 @@ int MkBuilder::filter_comb_cands(std::function filter) return n_removed; } -void MkBuilder::select_best_comb_cands(bool clear_m_tracks) +void MkBuilder::select_best_comb_cands(bool clear_m_tracks, bool remove_missing_hits) { if (clear_m_tracks) m_tracks.clear(); - export_best_comb_cands(m_tracks); + export_best_comb_cands(m_tracks, remove_missing_hits); } -void MkBuilder::export_best_comb_cands(TrackVec &out_vec) +void MkBuilder::export_best_comb_cands(TrackVec &out_vec, bool remove_missing_hits) { const EventOfCombCandidates &eoccs = m_event_of_comb_cands; out_vec.reserve(out_vec.size() + eoccs.m_size); @@ -323,7 +323,7 @@ void MkBuilder::export_best_comb_cands(TrackVec &out_vec) if ( ! eoccs[i].empty()) { const TrackCand &bcand = eoccs[i].front(); - out_vec.emplace_back( bcand.exportTrack() ); + out_vec.emplace_back( bcand.exportTrack(remove_missing_hits) ); } } } @@ -1961,7 +1961,10 @@ void MkBuilder::find_tracks_in_layers(CandCloner &cloner, MkFinder *mkfndr, // } // copy_out the propagated track params, errors only. - mkfndr->CopyOutParErr(eoccs.m_candidates, end - itrack, true); + // Do not, keep cands at last valid hit until actual update, + // this requires change to propagation flags used in MkFinder::UpdateWithLastHit() + // from intra-layer to inter-layer. + // mkfndr->CopyOutParErr(eoccs.m_candidates, end - itrack, true); dprint("make new candidates"); cloner.begin_iteration(); diff --git a/mkFit/MkBuilder.h b/mkFit/MkBuilder.h index 1e23d95c..3b6b5cad 100644 --- a/mkFit/MkBuilder.h +++ b/mkFit/MkBuilder.h @@ -160,8 +160,8 @@ class MkBuilder // XXX filter for rearranging cands that will / will not do backward search. - void select_best_comb_cands(bool clear_m_tracks=false); - void export_best_comb_cands(TrackVec &out_vec); + void select_best_comb_cands(bool clear_m_tracks=false, bool remove_missing_hits=false); + void export_best_comb_cands(TrackVec &out_vec, bool remove_missing_hits=false); void export_tracks(TrackVec &out_vec); void CompactifyHitStorageForBestCand(bool remove_seed_hits, int backward_fit_min_hits) diff --git a/mkFit/MkFinder.cc b/mkFit/MkFinder.cc index 446fac8e..47852487 100644 --- a/mkFit/MkFinder.cc +++ b/mkFit/MkFinder.cc @@ -1249,8 +1249,10 @@ void MkFinder::UpdateWithLastHit(const LayerOfHits &layer_of_hits, int N_proc, msPar.CopyIn(i, hit.posArray()); } + // See comment in MkBuilder::find_tracks_in_layer() about intra / inter flags used here + // for propagation to the hit. (*fnd_foos.m_update_param_foo)(Err[iP], Par[iP], Chg, msErr, msPar, - Err[iC], Par[iC], N_proc, Config::finding_intra_layer_pflags); + Err[iC], Par[iC], N_proc, Config::finding_inter_layer_pflags); } @@ -1527,6 +1529,7 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits, // Then we could avoid checking which layers actually do have hits. MPlexQF tmp_chi2; + MPlexQI no_mat_effs; float tmp_err[6] = { 666, 0, 666, 0, 0, 666 }; float tmp_pos[3]; @@ -1542,7 +1545,9 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits, const Hit *last_hit_ptr[NN]; #endif - int count = 0; + no_mat_effs.SetVal(0); + int done_count = 0; + int here_count = 0; for (int i = 0; i < N_proc; ++i) { while (CurNode[i] >= 0 && HoTNodeArr[ i ][ CurNode[i] ].m_hot.index < 0) @@ -1550,6 +1555,8 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits, CurNode[i] = HoTNodeArr[ i ][ CurNode[i] ].m_prev_idx; } + if (CurNode[i] < 0) ++done_count; + if (CurNode[i] >= 0 && HoTNodeArr[ i ][ CurNode[i] ].m_hot.layer == layer) { // Skip the overlap hits -- if they exist. @@ -1569,7 +1576,7 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits, #endif msErr.CopyIn(i, hit.errArray()); msPar.CopyIn(i, hit.posArray()); - ++count; + ++here_count; CurNode[i] = HoTNodeArr[ i ][ CurNode[i] ].m_prev_idx; } @@ -1579,6 +1586,7 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits, #ifdef DEBUG_BACKWARD_FIT last_hit_ptr[i] = nullptr; #endif + no_mat_effs[i] = 1; tmp_pos[0] = Par[iC](i, 0, 0); tmp_pos[1] = Par[iC](i, 1, 0); tmp_pos[2] = Par[iC](i, 2, 0); @@ -1587,20 +1595,21 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits, } } - if (count == 0) continue; + if (done_count == N_proc) break; + if (here_count == 0) continue; // ZZZ Could add missing hits here, only if there are any actual matches. if (LI.is_barrel()) { - PropagateTracksToHitR(msPar, N_proc, Config::backward_fit_pflags); + PropagateTracksToHitR(msPar, N_proc, Config::backward_fit_pflags, &no_mat_effs); kalmanOperation(KFO_Calculate_Chi2 | KFO_Update_Params, Err[iP], Par[iP], msErr, msPar, Err[iC], Par[iC], tmp_chi2, N_proc); } else { - PropagateTracksToHitZ(msPar, N_proc, Config::backward_fit_pflags); + PropagateTracksToHitZ(msPar, N_proc, Config::backward_fit_pflags, &no_mat_effs); kalmanOperationEndcap(KFO_Calculate_Chi2 | KFO_Update_Params, Err[iP], Par[iP], msErr, msPar, Err[iC], Par[iC], tmp_chi2, N_proc); diff --git a/mkFit/PropagationMPlex.cc b/mkFit/PropagationMPlex.cc index 059257c6..b90e9209 100644 --- a/mkFit/PropagationMPlex.cc +++ b/mkFit/PropagationMPlex.cc @@ -407,8 +407,9 @@ void helixAtRFromIterativeCCS(const MPlexLV& inPar, const MPlexQI& inChg, co void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const MPlexQI &inChg, const MPlexQF& msRad, - MPlexLS &outErr, MPlexLV& outPar, - const int N_proc, const PropagationFlags pflags) + MPlexLS &outErr, MPlexLV& outPar, + const int N_proc, const PropagationFlags pflags, + const MPlexQI *noMatEffPtr) { // bool debug = true; @@ -448,14 +449,17 @@ void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, MPlexQF hitsXi; MPlexQF propSign; #pragma omp simd - for (int n = 0; n < NN; ++n) + for (int n = 0; n < N_proc; ++n) { - const int zbin = getZbinME(outPar(n, 2, 0)); - const int rbin = getRbinME(msRad (n, 0, 0)); - - hitsRl(n, 0, 0) = (zbin>=0 && zbin=0 && rbin=0 && zbin=0 && rbinConstAt(n, 0, 0))) { + hitsRl(n, 0, 0) = 0.f; + hitsXi(n, 0, 0) = 0.f; + } else { + const int zbin = getZbinME(outPar(n, 2, 0)); + const int rbin = getRbinME(msRad (n, 0, 0)); + hitsRl(n, 0, 0) = (zbin>=0 && zbin=0 && rbin=0 && zbin=0 && rbinr0 ? 1. : -1.); @@ -505,8 +509,9 @@ void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const MPlexQI &inChg, const MPlexQF& msZ, - MPlexLS &outErr, MPlexLV& outPar, - const int N_proc, const PropagationFlags pflags) + MPlexLS &outErr, MPlexLV& outPar, + const int N_proc, const PropagationFlags pflags, + const MPlexQI *noMatEffPtr) { // debug = true; @@ -541,14 +546,17 @@ void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, MPlexQF hitsXi; MPlexQF propSign; #pragma omp simd - for (int n = 0; n < NN; ++n) + for (int n = 0; n < N_proc; ++n) { - const int zbin = getZbinME(msZ(n, 0, 0)); - const int rbin = getRbinME(std::hypot(outPar(n, 0, 0), outPar(n, 1, 0))); - - hitsRl(n, 0, 0) = (zbin>=0 && zbin=0 && rbin=0 && zbin=0 && rbinConstAt(n, 0, 0)) { + hitsRl(n, 0, 0) = 0.f; + hitsXi(n, 0, 0) = 0.f; + } else { + const int zbin = getZbinME(msZ(n, 0, 0)); + const int rbin = getRbinME(std::hypot(outPar(n, 0, 0), outPar(n, 1, 0))); + hitsRl(n, 0, 0) = (zbin>=0 && zbin=0 && rbin=0 && zbin=0 && rbinstd::abs(zin) ? 1. : -1.); diff --git a/mkFit/PropagationMPlex.h b/mkFit/PropagationMPlex.h index 8ac06bfd..a6873339 100644 --- a/mkFit/PropagationMPlex.h +++ b/mkFit/PropagationMPlex.h @@ -30,7 +30,8 @@ void propagateLineToRMPlex(const MPlexLS &psErr, const MPlexLV& psPar, void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const MPlexQI &inChg, const MPlexQF& msRad, MPlexLS &outErr, MPlexLV& outPar, - const int N_proc, const PropagationFlags pflags); + const int N_proc, const PropagationFlags pflags, + const MPlexQI *noMatEffPtr=nullptr); void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msRad, MPlexLV& outPar, MPlexLL& errorProp, @@ -44,7 +45,8 @@ void helixAtRFromIterativeCCS(const MPlexLV& inPar, const MPlexQI& inChg, const void propagateHelixToZMPlex(const MPlexLS &inErr, const MPlexLV& inPar, const MPlexQI &inChg, const MPlexQF& msZ, MPlexLS &outErr, MPlexLV& outPar, - const int N_proc, const PropagationFlags pflags); + const int N_proc, const PropagationFlags pflags, + const MPlexQI *noMatEffPtr=nullptr); void helixAtZ(const MPlexLV& inPar, const MPlexQI& inChg, const MPlexQF &msZ, MPlexLV& outPar, MPlexLL& errorProp, diff --git a/mkFit/buildtestMPlex.cc b/mkFit/buildtestMPlex.cc index 65ed228d..1abcbd65 100644 --- a/mkFit/buildtestMPlex.cc +++ b/mkFit/buildtestMPlex.cc @@ -709,7 +709,7 @@ void run_OneIteration(const TrackerInfo& trackerInfo, const IterationConfig &itc } } - builder.export_best_comb_cands(out_tracks); + builder.export_best_comb_cands(out_tracks, true); if (do_remove_duplicates) {