diff --git a/Config.cc b/Config.cc index cb266b66..c4e19471 100644 --- a/Config.cc +++ b/Config.cc @@ -105,6 +105,7 @@ namespace Config bool kludgeCmsHitErrors = false; bool backwardFit = false; + bool backwardSearch = true; bool includePCA = false; bool json_dump_before = false; diff --git a/Config.h b/Config.h index f7f709dd..22de32da 100644 --- a/Config.h +++ b/Config.h @@ -389,6 +389,7 @@ namespace Config extern bool kludgeCmsHitErrors; extern bool backwardFit; + extern bool backwardSearch; extern bool includePCA; // NAN and silly track parameter tracking options diff --git a/Geoms/CMS-2017.cc b/Geoms/CMS-2017.cc index 62481473..d57fb660 100644 --- a/Geoms/CMS-2017.cc +++ b/Geoms/CMS-2017.cc @@ -36,7 +36,7 @@ namespace sp.append_plan(47, false); sp.fill_plan(48, 53); // TID, 6 layers sp.fill_plan(54, 71); // TEC, 18 layers - sp.finalize_plan(); + sp.set_iterator_limits(2, 0); } { SteeringParams &sp = ic.m_steering_params[TrackerInfo::Reg_Transition_Neg]; @@ -51,7 +51,7 @@ namespace sp.fill_plan(48, 53); // TID, 6 layers sp.fill_plan(10, 17); // TOB, 8 layers sp.fill_plan(54, 71); // TEC, 18 layers - sp.finalize_plan(); + sp.set_iterator_limits(2, 0); } { SteeringParams &sp = ic.m_steering_params[TrackerInfo::Reg_Barrel]; @@ -59,9 +59,9 @@ namespace sp.fill_plan(0, 1, false, true); // bk-fit only sp.append_plan(2, true); // pickup-only sp.append_plan(3, false); - sp.fill_plan(4, 9); // TIB, 6 layers - sp.fill_plan(10, 17); // TOB, 8 layers - sp.finalize_plan(); + sp.fill_plan(4, 9); // TIB, 6 layers [ 4, 9] + sp.fill_plan(10, 17); // TOB, 8 layers [10, 17] + sp.set_iterator_limits(2, 0); } { SteeringParams &sp = ic.m_steering_params[TrackerInfo::Reg_Transition_Pos]; @@ -72,11 +72,11 @@ namespace sp.append_plan(18, false); sp.append_plan(19, false); sp.append_plan(20, false); - sp.fill_plan(4, 9); // TIB, 6 layers - sp.fill_plan(21, 26); // TID, 6 layers - sp.fill_plan(10, 17); // TOB, 8 layers - sp.fill_plan(27, 44); // TEC, 18 layers - sp.finalize_plan(); + sp.fill_plan(4, 9); // TIB, 6 layers [ 7, 12] + sp.fill_plan(21, 26); // TID, 6 layers [13, 18] + sp.fill_plan(10, 17); // TOB, 8 layers [19, 26] + sp.fill_plan(27, 44); // TEC, 18 layers [27, 44] + sp.set_iterator_limits(2, 0); } { SteeringParams &sp = ic.m_steering_params[TrackerInfo::Reg_Endcap_Pos]; @@ -86,12 +86,42 @@ namespace sp.append_plan(18, false); sp.append_plan(19, false); sp.append_plan(20, false); - sp.fill_plan(21, 26); // TID, 6 layers - sp.fill_plan(27, 44); // TEC, 18 layers - sp.finalize_plan(); + sp.fill_plan(21, 26); // TID, 6 layers [ 6, 11] + sp.fill_plan(27, 44); // TEC, 18 layers [12, 29] + sp.set_iterator_limits(2, 0); } } + void OverrideSteeringParams_Iter7(IterationConfig& ic) + { + ic.m_backward_search = true; + ic.m_backward_params = ic.m_params; + ic.m_backward_params.maxHolesPerCand = 2; + ic.m_backward_params.maxConsecHoles = 2; + // Remove pixel layers from FwdSearch, add them to BkwSearch + auto &spv = ic.m_steering_params; + spv[TrackerInfo::Reg_Endcap_Neg] .set_iterator_limits(8, 6, 19); + spv[TrackerInfo::Reg_Transition_Neg].set_iterator_limits(9, 7, 34); + spv[TrackerInfo::Reg_Barrel] .set_iterator_limits(6, 4, 8); + spv[TrackerInfo::Reg_Transition_Pos].set_iterator_limits(9, 7, 34); + spv[TrackerInfo::Reg_Endcap_Pos] .set_iterator_limits(8, 6, 19); + } + + void OverrideSteeringParams_Iter8(IterationConfig& ic) + { + ic.m_backward_search = true; + ic.m_backward_params = ic.m_params; + ic.m_backward_params.maxHolesPerCand = 2; + ic.m_backward_params.maxConsecHoles = 2; + // Remove pixel/tib/tid layers from FwdSearch, add them to BkwSearch/ + auto &spv = ic.m_steering_params; + spv[TrackerInfo::Reg_Endcap_Neg] .set_iterator_limits(12, 12, 24); + spv[TrackerInfo::Reg_Transition_Neg].set_iterator_limits(27, 19, 39); + spv[TrackerInfo::Reg_Barrel] .set_iterator_limits(12, 10, 14); + spv[TrackerInfo::Reg_Transition_Pos].set_iterator_limits(27, 19, 39); + spv[TrackerInfo::Reg_Endcap_Pos] .set_iterator_limits(12, 12, 24); + } + void SetupIterationParams(IterationParams& ip, unsigned int it=0) { if (it == 0) @@ -391,6 +421,7 @@ namespace fill_hit_selection_windows_params(ii[6]); ii[7].Clone(ii[0]); + OverrideSteeringParams_Iter7(ii[7]); SetupIterationParams(ii[7].m_params, 7); ii[7].set_iteration_index_and_track_algorithm(7, (int) TrackBase::TrackAlgorithm::pixelLessStep); ii[7].set_seed_cleaning_params(2.0, 0.135, 0.135, 0.135, 0.135, 0.135, 0.135, 0.135, 0.135); @@ -399,6 +430,7 @@ namespace fill_hit_selection_windows_params(ii[7]); ii[8].Clone(ii[0]); + OverrideSteeringParams_Iter8(ii[8]); SetupIterationParams(ii[8].m_params, 8); ii[8].set_iteration_index_and_track_algorithm(8, (int) TrackBase::TrackAlgorithm::tobTecStep); ii[8].set_seed_cleaning_params(2.0, 0.135, 0.135, 0.135, 0.135, 0.135, 0.135, 0.135, 0.135); diff --git a/mkFit/HitStructures.cc b/mkFit/HitStructures.cc index 305e7fcb..55ff1348 100644 --- a/mkFit/HitStructures.cc +++ b/mkFit/HitStructures.cc @@ -480,11 +480,11 @@ void CombCandidate::ImportSeed(const Track& seed, int region) { emplace_back(TrackCand(seed, this)); - m_state = CombCandidate::Dormant; - m_last_seed_layer = seed.getLastHitLyr(); + m_state = CombCandidate::Dormant; + m_pickup_layer = seed.getLastHitLyr(); #ifdef DUMPHITWINDOW - m_seed_algo = seed.algoint(); - m_seed_label = seed.label(); + m_seed_algo = seed.algoint(); + m_seed_label = seed.label(); #endif TrackCand &cand = back(); @@ -545,4 +545,64 @@ void CombCandidate::MergeCandsAndBestShortOne(const IterationParams& params, boo // assert(capacity() == (size_t)Config::maxCandsPerSeed); } +void CombCandidate::BeginBkwSearch() +{ + // The best candidate is assumed to be in position 0, after + // MergeCandsAndBestShortOne has been called. + // Other cands are dropped, their hits are not but they could be. + // This is to be called before backward-search to start with a single + // input candidate for backward cominatorial search. + // + // m_state and m_pickup_layer are also set. + + assert( ! empty() ); + + resize(1); + + TrackCand &tc = (*this)[0]; + + // Could compatify m_hots, removing all not used by element 0. + // Not trivial, as I have to start from back to get to the front + // and then reorder. Or have some semi-complicated scheme to + // put them at the back and then grow new hits in the front. + + m_state = Dormant; + m_pickup_layer = m_hots[0].m_hot.layer; + m_lastHitIdx_before_bkwsearch = tc.lastCcIndex(); + m_nInsideMinusOneHits_before_bkwsearch = tc.nInsideMinusOneHits(); + m_nTailMinusOneHits_before_bkwsearch = tc.nTailMinusOneHits(); + tc.setLastCcIndex(0); + tc.setNInsideMinusOneHits(0); + tc.setNTailMinusOneHits(0); +} + +void CombCandidate::EndBkwSearch() +{ + // MergeCandsAndBestShortOne() has already been called (from MkBuilder::FindXxx()). + // We have to fixup the best candidate. + + TrackCand &tc = (*this)[0]; + + int curr_idx = tc.lastCcIndex(); + if (curr_idx != 0) + { + int last_idx = -1, prev_idx; + do { + prev_idx = m_hots[curr_idx].m_prev_idx; + + m_hots[curr_idx].m_prev_idx = last_idx; + + last_idx = curr_idx; + curr_idx = prev_idx; + } while (prev_idx != -1); + } + + tc.setLastCcIndex(m_lastHitIdx_before_bkwsearch); + tc.setNInsideMinusOneHits(m_nInsideMinusOneHits_before_bkwsearch + tc.nInsideMinusOneHits()); + tc.setNTailMinusOneHits(m_nTailMinusOneHits_before_bkwsearch + tc.nTailMinusOneHits()); + m_lastHitIdx_before_bkwsearch = -1; + m_nInsideMinusOneHits_before_bkwsearch = -1; + m_nTailMinusOneHits_before_bkwsearch = -1; +} + } // end namespace mkfit diff --git a/mkFit/HitStructures.h b/mkFit/HitStructures.h index 156848d7..e807c672 100644 --- a/mkFit/HitStructures.h +++ b/mkFit/HitStructures.h @@ -477,7 +477,7 @@ inline float getScoreCand(const TrackCand& cand1, bool penalizeTailMissHits=fals int nfoundhits = cand1.nFoundHits(); int noverlaphits = cand1.nOverlapHits(); int nmisshits = cand1.nInsideMinusOneHits(); - float ntailmisshits = penalizeTailMissHits ? cand1.nTailMinusOneHits() : 0; + int ntailmisshits = penalizeTailMissHits ? cand1.nTailMinusOneHits() : 0; float pt = cand1.pT(); float chi2 = cand1.chi2(); // Do not allow for chi2<0 in score calculation @@ -494,8 +494,12 @@ class CombCandidate : public std::vector enum SeedState_e { Dormant = 0, Finding, Finished }; TrackCand m_best_short_cand; - SeedState_e m_state = Dormant; - int m_last_seed_layer = -1; + SeedState_e m_state : 8; + int m_pickup_layer : 16; + short int m_lastHitIdx_before_bkwsearch = -1; + short int m_nInsideMinusOneHits_before_bkwsearch = -1; + short int m_nTailMinusOneHits_before_bkwsearch = -1; + #ifdef DUMPHITWINDOW int m_seed_algo = 0; int m_seed_label = 0; @@ -506,16 +510,19 @@ class CombCandidate : public std::vector std::vector m_overlap_hits; // XXXX HitMatchPair could be a member in TrackCand - CombCandidate() - { - } + CombCandidate() : + m_state(Dormant), m_pickup_layer(-1) + {} // Need this so resize of EventOfCombinedCandidates::m_candidates can reuse vectors used here. CombCandidate(CombCandidate&& o) : std::vector(std::move(o)), m_best_short_cand(std::move(o.m_best_short_cand)), m_state(o.m_state), - m_last_seed_layer(o.m_last_seed_layer), + m_pickup_layer(o.m_pickup_layer), + m_lastHitIdx_before_bkwsearch(o.m_lastHitIdx_before_bkwsearch), + m_nInsideMinusOneHits_before_bkwsearch(o.m_nInsideMinusOneHits_before_bkwsearch), + m_nTailMinusOneHits_before_bkwsearch(o.m_nTailMinusOneHits_before_bkwsearch), #ifdef DUMPHITWINDOW m_seed_algo(o.m_seed_algo), m_seed_label(o.m_seed_label), @@ -540,7 +547,10 @@ class CombCandidate : public std::vector std::vector::operator=( std::move(o) ); m_best_short_cand = std::move(o.m_best_short_cand); m_state = o.m_state; - m_last_seed_layer = o.m_last_seed_layer; + m_pickup_layer = o.m_pickup_layer; + m_lastHitIdx_before_bkwsearch = o.m_lastHitIdx_before_bkwsearch; + m_nInsideMinusOneHits_before_bkwsearch = o.m_nInsideMinusOneHits_before_bkwsearch; + m_nTailMinusOneHits_before_bkwsearch = o.m_nTailMinusOneHits_before_bkwsearch; #ifdef DUMPHITWINDOW m_seed_algo = o.m_seed_algo; m_seed_label = o.m_seed_label; @@ -561,6 +571,8 @@ class CombCandidate : public std::vector m_best_short_cand.setScore( getScoreWorstPossible() ); + // state and pickup_layer set in ImportSeed. + // expected_num_hots is different for CloneEngine and Std, especially as long as we // instantiate all candidates before purging them. // ce: N_layer * N_cands ~~ 20 * 6 = 120 @@ -591,6 +603,9 @@ class CombCandidate : public std::vector } void MergeCandsAndBestShortOne(const IterationParams¶ms, bool update_score, bool sort_cands); + + void BeginBkwSearch(); + void EndBkwSearch(); }; //============================================================================== @@ -689,6 +704,9 @@ class EventOfCombCandidates ++m_size; } + + void BeginBkwSearch() { for (int i=0; i */ m_region_order, // /* vector */ m_steering_params, diff --git a/mkFit/IterationConfig.h b/mkFit/IterationConfig.h index a403479b..427394d4 100644 --- a/mkFit/IterationConfig.h +++ b/mkFit/IterationConfig.h @@ -155,8 +155,11 @@ class IterationConfig bool m_require_quality_filter = false; bool m_require_dupclean_tight = false; + bool m_backward_search = false; + // Iteration parameters (could be a ptr) IterationParams m_params; + IterationParams m_backward_params; int m_n_regions; std::vector m_region_order; @@ -239,6 +242,7 @@ class IterationConfig m_n_regions = nreg; m_region_order.resize(nreg); m_steering_params.resize(nreg); + for (int i = 0; i < nreg; ++i) m_steering_params[i].m_region = i; m_layer_configs.resize(nlay); } diff --git a/mkFit/MkBuilder.cc b/mkFit/MkBuilder.cc index a3d89468..ec91b623 100644 --- a/mkFit/MkBuilder.cc +++ b/mkFit/MkBuilder.cc @@ -303,8 +303,10 @@ int MkBuilder::filter_comb_cands(std::function filter) return n_removed; } -void MkBuilder::select_best_comb_cands() +void MkBuilder::select_best_comb_cands(bool clear_m_tracks) { + if (clear_m_tracks) + m_tracks.clear(); export_best_comb_cands(m_tracks); } @@ -902,7 +904,7 @@ void MkBuilder::quality_print() } } -void MkBuilder::track_print(Track &t, const char* pref) +void MkBuilder::track_print(const Track &t, const char* pref) { printf("%s with q=%+i pT=%7.3f eta=% 7.3f nHits=%2d label=%4d\nState:\n", pref, t.charge(), t.pT(), t.momEta(), t.nFoundHits(), t.label()); @@ -1126,7 +1128,11 @@ void MkBuilder::seed_post_cleaning(TrackVec &tv, const bool fix_silly_seeds, con break; } } - if (tv.size() != 1) printf("Preselect seed with label %d - NOT FOUND. Running on full event.\n", SELECT_SEED_LABEL); + if (tv.size() != 1) + { + printf("Preselect seed with label %d - NOT FOUND. Cleaning out seeds.\n", SELECT_SEED_LABEL); + tv.clear(); + } } #endif @@ -1287,7 +1293,7 @@ void MkBuilder::find_tracks_load_seeds_BH(const TrackVec &in_seeds) } -void MkBuilder::FindTracksBestHit() +void MkBuilder::FindTracksBestHit(SteeringParams::IterationType_e iteration_dir) { // bool debug = true; @@ -1296,6 +1302,11 @@ void MkBuilder::FindTracksBestHit() tbb::parallel_for_each(m_job->regions_begin(), m_job->regions_end(), [&](int region) { + if (iteration_dir == SteeringParams::IT_BkwSearch && ! m_job->steering_params(region).has_bksearch_plan()) { + printf("No backward search plan for region %d\n", region); + return; + } + // XXXXXX Select endcap / barrel only ... // if (region != TrackerInfo::Reg_Endcap_Neg && region != TrackerInfo::Reg_Endcap_Pos) // if (region != TrackerInfo::Reg_Barrel) @@ -1332,21 +1343,24 @@ void MkBuilder::FindTracksBestHit() } int curr_tridx = 0; - auto layer_plan_it = st_par.finding_begin(); + auto layer_plan_it = st_par.make_iterator(iteration_dir); - assert( layer_plan_it->m_pickup_only ); + dprintf("Made iterator for %d, first layer=%d ... end layer=%d\n", + iteration_dir, layer_plan_it.layer(), layer_plan_it.last_layer()); - int curr_layer = layer_plan_it->m_layer; + assert( layer_plan_it.is_pickup_only() ); + + int curr_layer = layer_plan_it.layer(); mkfndr->Stopped.SetVal(0); // Loop over layers, starting from after the seed. // Consider inverting loop order and make layer outer, need to // trade off hit prefetching with copy-out of candidates. - while (++layer_plan_it != st_par.finding_end()) + while (++layer_plan_it) { prev_layer = curr_layer; - curr_layer = layer_plan_it->m_layer; + curr_layer = layer_plan_it.layer(); mkfndr->Setup(m_job->m_iter_config.m_params, m_job->m_iter_config.m_layer_configs[curr_layer], m_job->get_mask_for_layer(curr_layer)); @@ -1372,7 +1386,7 @@ void MkBuilder::FindTracksBestHit() } } - if (layer_plan_it->m_pickup_only) continue; + if (layer_plan_it.is_pickup_only()) continue; dcall(pre_prop_print(curr_layer, mkfndr.get())); @@ -1463,7 +1477,7 @@ int MkBuilder::find_tracks_unroll_candidates(std::vector> & s { CombCandidate &ccand = m_event_of_comb_cands[iseed]; - if (ccand.m_state == CombCandidate::Dormant && ccand.m_last_seed_layer == prev_layer) + if (ccand.m_state == CombCandidate::Dormant && ccand.m_pickup_layer == prev_layer) { ccand.m_state = CombCandidate::Finding; } @@ -1577,7 +1591,7 @@ void MkBuilder::find_tracks_handle_missed_layers(MkFinder *mkfndr, const LayerIn // FindTracksCombinatorial: Standard TBB //------------------------------------------------------------------------------ -void MkBuilder::FindTracksStandard() +void MkBuilder::FindTracksStandard(SteeringParams::IterationType_e iteration_dir) { // debug = true; @@ -1586,6 +1600,11 @@ void MkBuilder::FindTracksStandard() tbb::parallel_for_each(m_job->regions_begin(), m_job->regions_end(), [&](int region) { + if (iteration_dir == SteeringParams::IT_BkwSearch && ! m_job->steering_params(region).has_bksearch_plan()) { + printf("No backward search plan for region %d\n", region); + return; + } + const TrackerInfo &trk_info = m_job->m_trk_info; const SteeringParams &st_par = m_job->steering_params(region); const IterationParams ¶ms = m_job->params(); @@ -1615,21 +1634,28 @@ void MkBuilder::FindTracksStandard() std::vector> seed_cand_idx; seed_cand_idx.reserve(n_seeds * params.maxCandsPerSeed); - auto layer_plan_it = st_par.finding_begin(); + auto layer_plan_it = st_par.make_iterator(iteration_dir); - assert( layer_plan_it->m_pickup_only ); + dprintf("Made iterator for %d, first layer=%d ... end layer=%d\n", + iteration_dir, layer_plan_it.layer(), layer_plan_it.last_layer()); - int curr_layer = layer_plan_it->m_layer, prev_layer; + assert( layer_plan_it.is_pickup_only() ); + + int curr_layer = layer_plan_it.layer(), prev_layer; dprintf("\nMkBuilder::FindTracksStandard region=%d, seed_pickup_layer=%d, first_layer=%d\n", - region, curr_layer, (layer_plan_it + 1)->m_layer); + region, curr_layer, layer_plan_it.next_layer()); + + auto &iter_params = (iteration_dir == SteeringParams::IT_BkwSearch) ? + m_job->m_iter_config.m_backward_params : + m_job->m_iter_config.m_params; // Loop over layers, starting from after the seed. - while (++layer_plan_it != st_par.finding_end()) + while (++layer_plan_it) { prev_layer = curr_layer; - curr_layer = layer_plan_it->m_layer; - mkfndr->Setup(m_job->m_iter_config.m_params, m_job->m_iter_config.m_layer_configs[curr_layer], + curr_layer = layer_plan_it.layer(); + mkfndr->Setup(iter_params, m_job->m_iter_config.m_layer_configs[curr_layer], m_job->get_mask_for_layer(curr_layer)); dprintf("\n* Processing layer %d\n", curr_layer); @@ -1639,9 +1665,9 @@ void MkBuilder::FindTracksStandard() const FindingFoos &fnd_foos = FindingFoos::get_finding_foos(layer_info.is_barrel()); int theEndCand = find_tracks_unroll_candidates(seed_cand_idx, start_seed, end_seed, - prev_layer, layer_plan_it->m_pickup_only); + prev_layer, layer_plan_it.is_pickup_only()); - if (layer_plan_it->m_pickup_only || theEndCand == 0) continue; + if (layer_plan_it.is_pickup_only() || theEndCand == 0) continue; // vectorized loop for (int itrack = 0; itrack < theEndCand; itrack += NN) @@ -1768,7 +1794,7 @@ void MkBuilder::FindTracksStandard() // FindTracksCombinatorial: CloneEngine TBB //------------------------------------------------------------------------------ -void MkBuilder::FindTracksCloneEngine() +void MkBuilder::FindTracksCloneEngine(SteeringParams::IterationType_e iteration_dir) { // debug = true; @@ -1777,6 +1803,11 @@ void MkBuilder::FindTracksCloneEngine() tbb::parallel_for_each(m_job->regions_begin(), m_job->regions_end(), [&](int region) { + if (iteration_dir == SteeringParams::IT_BkwSearch && ! m_job->steering_params(region).has_bksearch_plan()) { + printf("No backward search plan for region %d\n", region); + return; + } + const RegionOfSeedIndices rosi(m_seedEtaSeparators, region); // adaptive seeds per task based on the total estimated amount of work to divide among all threads @@ -1792,7 +1823,7 @@ void MkBuilder::FindTracksCloneEngine() cloner->Setup(m_job->params()); // loop over layers - find_tracks_in_layers(*cloner, mkfndr.get(), seeds.begin(), seeds.end(), region); + find_tracks_in_layers(*cloner, mkfndr.get(), iteration_dir, seeds.begin(), seeds.end(), region); cloner->Release(); }); @@ -1802,6 +1833,7 @@ void MkBuilder::FindTracksCloneEngine() } void MkBuilder::find_tracks_in_layers(CandCloner &cloner, MkFinder *mkfndr, + SteeringParams::IterationType_e iteration_dir, const int start_seed, const int end_seed, const int region) { EventOfCombCandidates &eoccs = m_event_of_comb_cands; @@ -1824,24 +1856,31 @@ void MkBuilder::find_tracks_in_layers(CandCloner &cloner, MkFinder *mkfndr, // Note that we do a final pass with curr_layer = -1 to update parameters // and output final tracks. - auto layer_plan_it = st_par.finding_begin(); + auto layer_plan_it = st_par.make_iterator(iteration_dir); + + dprintf("Made iterator for %d, first layer=%d ... end layer=%d\n", + iteration_dir, layer_plan_it.layer(), layer_plan_it.last_layer()); - assert( layer_plan_it->m_pickup_only ); + assert( layer_plan_it.is_pickup_only() ); - int curr_layer = layer_plan_it->m_layer, prev_layer; + int curr_layer = layer_plan_it.layer(), prev_layer; dprintf("\nMkBuilder::find_tracks_in_layers region=%d, seed_pickup_layer=%d, first_layer=%d; start_seed=%d, end_seed=%d\n", - region, curr_layer, (layer_plan_it + 1)->m_layer, start_seed, end_seed); + region, curr_layer, layer_plan_it.next_layer(), start_seed, end_seed); + + auto &iter_params = (iteration_dir == SteeringParams::IT_BkwSearch) ? + m_job->m_iter_config.m_backward_params : + m_job->m_iter_config.m_params; // Loop over layers according to plan. - while (++layer_plan_it != st_par.finding_end()) + while (++layer_plan_it) { prev_layer = curr_layer; - curr_layer = layer_plan_it->m_layer; - mkfndr->Setup(m_job->m_iter_config.m_params, m_job->m_iter_config.m_layer_configs[curr_layer], + curr_layer = layer_plan_it.layer(); + mkfndr->Setup(iter_params, m_job->m_iter_config.m_layer_configs[curr_layer], m_job->get_mask_for_layer(curr_layer)); - const bool pickup_only = layer_plan_it->m_pickup_only; + const bool pickup_only = layer_plan_it.is_pickup_only(); dprintf("\n\n* Processing layer %d, %s\n\n", curr_layer, pickup_only ? "pickup only" : "full finding"); @@ -2093,6 +2132,9 @@ void MkBuilder::fit_cands(MkFinder *mkfndr, int start_cand, int end_cand, int re // Check if we need to fragment this for SlurpIn to work. // Would actually prefer to do memory allocator for HoTNode storage. + // This has been "fixed" by copying Cands back into original container, not swapping the contents + // with vectors created in another thread (and thus not in the same memory pool, apparently), see + // CandCloner::ProcessSeedRange(). Standard building does not have this problem. /* step = NN; { diff --git a/mkFit/MkBuilder.h b/mkFit/MkBuilder.h index 6f8092b9..5d602a6a 100644 --- a/mkFit/MkBuilder.h +++ b/mkFit/MkBuilder.h @@ -153,10 +153,13 @@ class MkBuilder void import_seeds(const TrackVec &in_seeds, std::function insert_seed); int filter_comb_cands(std::function filter); - void select_best_comb_cands(); + void select_best_comb_cands(bool clear_m_tracks=false); void export_best_comb_cands(TrackVec &out_vec); void export_tracks(TrackVec &out_vec); + void BeginBkwSearch() { m_event_of_comb_cands.BeginBkwSearch(); } + void EndBkwSearch() { m_event_of_comb_cands.EndBkwSearch(); } + // MIMI hack to export tracks for BH const TrackVec& ref_tracks() const { return m_tracks; } TrackVec& ref_tracks_nc() { return m_tracks; } @@ -171,7 +174,7 @@ class MkBuilder void quality_reset(); void quality_process(Track& tkcand, const int itrack, std::map & cmsswLabelToPos); void quality_print(); - void track_print(Track &t, const char* pref); + void track_print(const Track &t, const char* pref); void quality_store_tracks(TrackVec & tracks); @@ -201,6 +204,7 @@ class MkBuilder const int itrack, const int end); void find_tracks_in_layers(CandCloner &cloner, MkFinder *mkfndr, + SteeringParams::IterationType_e iteration_dir, const int start_seed, const int end_seed, const int region); // -------- @@ -208,9 +212,9 @@ class MkBuilder void PrepareSeeds(); - void FindTracksBestHit(); - void FindTracksStandard(); - void FindTracksCloneEngine(); + void FindTracksBestHit(SteeringParams::IterationType_e iteration_dir=SteeringParams::IT_FwdSearch); + void FindTracksStandard(SteeringParams::IterationType_e iteration_dir=SteeringParams::IT_FwdSearch); + void FindTracksCloneEngine(SteeringParams::IterationType_e iteration_dir=SteeringParams::IT_FwdSearch); void BackwardFitBH(); void fit_cands_BH(MkFinder *mkfndr, int start_cand, int end_cand, int region); diff --git a/mkFit/MkFinder.cc b/mkFit/MkFinder.cc index 61f41253..f43223f6 100644 --- a/mkFit/MkFinder.cc +++ b/mkFit/MkFinder.cc @@ -1372,9 +1372,9 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits, float tmp_err[6] = { 666, 0, 666, 0, 0, 666 }; float tmp_pos[3]; - for (auto lp_iter = st_par.m_layer_plan.rbegin(); lp_iter != st_par.m_layer_plan.rend(); ++lp_iter) + for (auto lp_iter = st_par.make_iterator(SteeringParams::IT_BkwFit); lp_iter.is_valid(); ++lp_iter) { - const int layer = lp_iter->m_layer; + const int layer = lp_iter.layer(); const LayerOfHits &L = eventofhits.m_layers_of_hits[layer]; const LayerInfo &LI = * L.m_layer_info; @@ -1444,6 +1444,16 @@ void MkFinder::BkFitFitTracks(const EventOfHits & eventofhits, Err[iP], Par[iP], msErr, msPar, Err[iC], Par[iC], tmp_chi2, N_proc); } + //fixup invpt sign and charge + for (int n = 0; n < N_proc; ++n) + { + if (Par[iC].At(n,3,0) < 0) + { + Chg.At(n, 0, 0) = -Chg.At(n, 0, 0); + Par[iC].At(n,3,0) = -Par[iC].At(n,3,0); + } + } + for (int i = 0; i < N_proc; ++i) { #ifdef DEBUG_BACKWARD_FIT diff --git a/mkFit/SteeringParams.h b/mkFit/SteeringParams.h index 31e80c60..8a91a8fd 100644 --- a/mkFit/SteeringParams.h +++ b/mkFit/SteeringParams.h @@ -2,6 +2,7 @@ #define SteeringParams_h #include +#include namespace mkfit { @@ -18,13 +19,23 @@ struct LayerControl // int m_on_miss_jump_to = -999; // int m_on_hit_jump_to = -999; - bool m_pickup_only; // do not propagate to this layer and process hits, pickup seeds only. - bool m_backfit_only; // layer only used in backward fit. + // XXX To Be removed << FNORD + bool m_pickup_only; // do not propagate to this layer and process hits, pickup seeds only. + bool m_bkfit_only; // layer only used in backward fit. + bool m_bksearch_only; // layer only used in backward search. + + bool skip_on_fwd() const { return ! m_pickup_only; } + bool skip_on_bkfit() const { return m_bksearch_only; } + bool use_in_bkwsearch() const { return m_bksearch_only || m_bkfit_only || m_pickup_only; } + // FNORD //---------------------------------------------------------------------------- - LayerControl(int lay, bool pu_only=false, bool bf_only=false) : - m_layer(lay), m_pickup_only(pu_only), m_backfit_only(bf_only) {} + LayerControl(int lay, bool pu_only=false, bool bf_only=false, bool bs_only=false) : + m_layer(lay), m_pickup_only(pu_only), m_bkfit_only(bf_only), m_bksearch_only(bs_only) {} + + void fixup(bool pu_only, bool bf_only, bool bs_only) + { m_pickup_only = pu_only; m_bkfit_only = bf_only; m_bksearch_only = bs_only; } }; @@ -35,9 +46,73 @@ struct LayerControl class SteeringParams { public: - std::vector m_layer_plan; - std::vector::iterator m_begin_for_finding; - std::vector::iterator m_end_for_finding; + enum IterationType_e { IT_FwdSearch, IT_BkwFit, IT_BkwSearch }; + + class iterator + { + friend class SteeringParams; + + const SteeringParams &m_steering_params; + IterationType_e m_type; + int m_cur_index = -1; + int m_end_index = -1; + + iterator(const SteeringParams& sp, IterationType_e t) : + m_steering_params(sp), m_type(t) + {} + +public: + const LayerControl& layer_control() const { return m_steering_params.m_layer_plan[m_cur_index]; } + int layer() const { return layer_control().m_layer; } + int index() const { return m_cur_index; } + int region() const { return m_steering_params.m_region; } + + bool is_valid() const { return m_cur_index != -1; } + + const LayerControl& operator->() const { return layer_control(); } + + bool is_pickup_only() const + { + if (m_type == IT_FwdSearch) return m_cur_index == m_steering_params.m_fwd_search_pickup; + else if (m_type == IT_BkwSearch) return m_cur_index == m_steering_params.m_bkw_search_pickup; + else throw std::runtime_error("invalid iteration type"); + } + + bool operator++() + { + if ( ! is_valid()) return false; + if (m_type == IT_FwdSearch) + { + if (++m_cur_index == m_end_index) + m_cur_index = -1; + } + else + { + if (--m_cur_index == m_end_index) + m_cur_index = -1; + } + return is_valid(); + } + + // Functions for debug printouts + int end_index() const { return m_end_index; } + int next_layer() const { + if (m_type == IT_FwdSearch) return m_steering_params.m_layer_plan[m_cur_index + 1].m_layer; + else return m_steering_params.m_layer_plan[m_cur_index - 1].m_layer; + } + int last_layer() const { + if (m_type == IT_FwdSearch) return m_steering_params.m_layer_plan[m_end_index - 1].m_layer; + else return m_steering_params.m_layer_plan[m_end_index + 1].m_layer; + } + }; + + std::vector m_layer_plan; + + int m_region; + + int m_fwd_search_pickup = 0; + int m_bkw_fit_last = 0; + int m_bkw_search_pickup = -1; //---------------------------------------------------------------------------- @@ -55,18 +130,52 @@ class SteeringParams void fill_plan(int first, int last, bool pu_only=false, bool bf_only=false) { - for (int i = first; i <= last; ++i) append_plan(i, pu_only, bf_only); + for (int i = first; i <= last; ++i) append_plan(i, pu_only, bf_only); } - void finalize_plan() + void fixup_plan(int first_idx, int last_idx, bool pu_only, bool bf_only, bool bs_only) { - m_begin_for_finding = m_layer_plan.begin(); - while (m_begin_for_finding->m_backfit_only) ++m_begin_for_finding; - m_end_for_finding = m_layer_plan.end(); + for (int i = first_idx; i <= last_idx; ++i) + m_layer_plan[i].fixup(pu_only, bf_only, bs_only); } - std::vector::iterator finding_begin() const { return m_begin_for_finding; } - std::vector::iterator finding_end() const { return m_end_for_finding; } + void set_iterator_limits(int fwd_search_pu, int bkw_fit_last, int bkw_search_pu=-1) + { + m_fwd_search_pickup = fwd_search_pu; + m_bkw_fit_last = bkw_fit_last; + m_bkw_search_pickup = bkw_search_pu; + } + + bool has_bksearch_plan() const + { + return m_bkw_search_pickup != -1; + } + + iterator make_iterator(IterationType_e type) const + { + iterator it(*this, type); + + if (type == IT_FwdSearch) + { + it.m_cur_index = m_fwd_search_pickup; + it.m_end_index = m_layer_plan.size(); + } + else if (type == IT_BkwFit) + { + it.m_cur_index = m_layer_plan.size() - 1; + it.m_end_index = m_bkw_fit_last - 1; + } + else if (type == IT_BkwSearch) + { + it.m_cur_index = m_bkw_search_pickup; + it.m_end_index = -1; + } + else throw std::invalid_argument("unknown iteration type"); + + if ( ! it.is_valid()) throw std::runtime_error("invalid iterator constructed"); + + return it; + } }; } // end namespace mkfit diff --git a/mkFit/buildtestMPlex.cc b/mkFit/buildtestMPlex.cc index 3d118aae..a06f6df9 100644 --- a/mkFit/buildtestMPlex.cc +++ b/mkFit/buildtestMPlex.cc @@ -437,7 +437,7 @@ std::vector runBtpCe_MultiIter(Event& ev, const EventOfHits &eoh, MkBuil TrackVec seeds_used; TrackVec seeds1; - + unsigned int algorithms[]={ 4,22,23,5,24,7,8,9,10,6 }; //9 iterations if (validation_on) @@ -489,13 +489,13 @@ std::vector runBtpCe_MultiIter(Event& ev, const EventOfHits &eoh, MkBuil if ( ! itconf.m_require_quality_filter) StdSeq::clean_cms_seedtracks_iter(&seeds, itconf); - + + builder.seed_post_cleaning(seeds, true, true); + // Add protection in case no seeds are found for iteration if (seeds.size() <= 0) continue; - builder.seed_post_cleaning(seeds, true, true); - builder.find_tracks_load_seeds(seeds); double time = dtime(); @@ -532,17 +532,24 @@ std::vector runBtpCe_MultiIter(Event& ev, const EventOfHits &eoh, MkBuil if (Config::backwardFit) { // a) TrackVec version: - builder.BackwardFitBH(); + // builder.BackwardFitBH(); - StdSeq::find_and_remove_duplicates(builder.ref_tracks_nc(), itconf); + // b) Version that runs on CombCand / TrackCand + builder.BackwardFit(); - builder.export_tracks(ev.fitTracks_); + if (Config::backwardSearch && itconf.m_backward_search) + { + builder.BeginBkwSearch(); + builder.FindTracksCloneEngine(SteeringParams::IT_BkwSearch); + builder.EndBkwSearch(); + } - // b) Version that runs on CombCand / TrackCand - // builder.BackwardFit(); - // builder.quality_store_tracks(ev.fitTracks_); + builder.select_best_comb_cands(true); // true -> clear m_tracks as they were already filled once above + + StdSeq::find_and_remove_duplicates(builder.ref_tracks_nc(), itconf); + builder.export_tracks(ev.fitTracks_); } - + builder.end_event(); } @@ -639,19 +646,19 @@ void run_OneIteration(const TrackerInfo& trackerInfo, const IterationConfig &itc { return StdSeq::qfilter_n_hits_pixseed(t, 3); }); } - builder.select_best_comb_cands(); - if (do_backward_fit) { - // a) BackwardFitBH works on MkBuilder::m_tracks - builder.BackwardFitBH(); + builder.BackwardFit(); - // b) Version that runs on CombCand / TrackCand - // builder.BackwardFit(); - // builder.export_best_comb_cands(out_tracks); + if (itconf.m_backward_search) + { + builder.BeginBkwSearch(); + builder.FindTracksCloneEngine(SteeringParams::IT_BkwSearch); + builder.EndBkwSearch(); + } } - builder.export_tracks(out_tracks); + builder.export_best_comb_cands(out_tracks); if (do_remove_duplicates) { diff --git a/mkFit/mkFit.cc b/mkFit/mkFit.cc index a37fc1ad..2e407413 100644 --- a/mkFit/mkFit.cc +++ b/mkFit/mkFit.cc @@ -347,7 +347,7 @@ void test_standard() StdSeq::LoadDeads(eoh, deadvectors); } - double t_best[NT] = {0}, t_cur[NT]; + double t_best[NT] = {0}, t_cur[NT] = {0}; std::vector t_cur_iter; simtrackstot += ev.simTracks_.size(); seedstot += ev.seedTracks_.size(); @@ -568,6 +568,8 @@ int main(int argc, const char *argv[]) " --use-phiq-arr use phi-Q arrays in select hit indices (def: %s)\n" " --kludge-cms-hit-errors make sure err(xy) > 15 mum, err(z) > 30 mum (def: %s)\n" " --backward-fit perform backward fit during building (def: %s)\n" + " --no-backward-search do not do backward search after backward fit\n" + " (def: do search if backward-fit is enabled and available in given iteration)\n" " --include-pca do the backward fit to point of closest approach, does not imply '--backward-fit' (def: %s)\n" "\n----------------------------------------------------------------------------------------------------------\n\n" "Validation options\n\n" @@ -911,6 +913,10 @@ int main(int argc, const char *argv[]) { Config::backwardFit = true; } + else if(*i == "--no-backward-search") + { + Config::backwardSearch = false; + } else if(*i == "--include-pca") { Config::includePCA = true;