MkFitCore/src/MkBuilder.cc

0001 #include <memory>
0002 #include <limits>
0003 #include <algorithm>
0004
0005 #include "RecoTracker/MkFitCore/interface/cms_common_macros.h"
0006
0007 #include "RecoTracker/MkFitCore/interface/MkBuilder.h"
0008 #include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
0009 #include "RecoTracker/MkFitCore/interface/binnor.h"
0010
0011 #include "Pool.h"
0012 #include "CandCloner.h"
0013 #include "FindingFoos.h"
0014 #include "MkFitter.h"
0015 #include "MkFinder.h"
0016
0017 #ifdef MKFIT_STANDALONE
0018 #include "RecoTracker/MkFitCore/standalone/Event.h"
0019 #endif
0020
0021 //#define DEBUG
0022 #include "Debug.h"
0023 //#define DEBUG_FINAL_FIT
0024
0025 #include "oneapi/tbb/parallel_for.h"
0026 #include "oneapi/tbb/parallel_for_each.h"
0027
0028 namespace mkfit {
0029
0030   //==============================================================================
0031   // Execution context -- Pools of helper objects
0032   //==============================================================================
0033
0034   struct ExecutionContext {
0035     ExecutionContext() = default;
0036     ~ExecutionContext() = default;
0037
0038     Pool<CandCloner> m_cloners;
0039     Pool<MkFitter> m_fitters;
0040     Pool<MkFinder> m_finders;
0041
0042     void populate(int n_thr) {
0043       m_cloners.populate(n_thr - m_cloners.size());
0044       m_fitters.populate(n_thr - m_fitters.size());
0045       m_finders.populate(n_thr - m_finders.size());
0046     }
0047
0048     void clear() {
0049       m_cloners.clear();
0050       m_fitters.clear();
0051       m_finders.clear();
0052     }
0053   };
0054
0055   CMS_SA_ALLOW ExecutionContext g_exe_ctx;
0056
0057 }  // end namespace mkfit
0058
0059 //------------------------------------------------------------------------------
0060
0061 namespace {
0062   using namespace mkfit;
0063
0064   // Range of indices processed within one iteration of a TBB parallel_for.
0065   struct RangeOfSeedIndices {
0066     int m_rng_beg, m_rng_end;
0067     int m_beg, m_end;
0068
0069     RangeOfSeedIndices(int rb, int re) : m_rng_beg(rb), m_rng_end(re) { reset(); }
0070
0071     void reset() {
0072       m_end = m_rng_beg;
0073       next_chunk();
0074     }
0075
0076     bool valid() const { return m_beg < m_rng_end; }
0077
0078     int n_proc() const { return m_end - m_beg; }
0079
0080     void next_chunk() {
0081       m_beg = m_end;
0082       m_end = std::min(m_end + NN, m_rng_end);
0083     }
0084
0085     RangeOfSeedIndices &operator++() {
0086       next_chunk();
0087       return *this;
0088     }
0089   };
0090
0091   // Region of seed indices processed in a single TBB parallel for.
0092   struct RegionOfSeedIndices {
0093     int m_reg_beg, m_reg_end, m_vec_cnt;
0094
0095     RegionOfSeedIndices(const std::vector<int> &seedEtaSeparators, int region) {
0096       m_reg_beg = (region == 0) ? 0 : seedEtaSeparators[region - 1];
0097       m_reg_end = seedEtaSeparators[region];
0098       m_vec_cnt = (m_reg_end - m_reg_beg + NN - 1) / NN;
0099     }
0100
0101     int count() const { return m_reg_end - m_reg_beg; }
0102
0103     tbb::blocked_range<int> tbb_blk_rng_std(int thr_hint = -1) const {
0104       if (thr_hint < 0)
0105         thr_hint = Config::numSeedsPerTask;
0106       return tbb::blocked_range<int>(m_reg_beg, m_reg_end, thr_hint);
0107     }
0108
0109     tbb::blocked_range<int> tbb_blk_rng_vec() const {
0110       return tbb::blocked_range<int>(0, m_vec_cnt, std::max(1, Config::numSeedsPerTask / NN));
0111     }
0112
0113     RangeOfSeedIndices seed_rng(const tbb::blocked_range<int> &i) const {
0114       return RangeOfSeedIndices(m_reg_beg + NN * i.begin(), std::min(m_reg_beg + NN * i.end(), m_reg_end));
0115     }
0116   };
0117
0118 #ifdef DEBUG
0119   void pre_prop_print(int ilay, MkBase *fir) {
0120     const float pt = 1.f / fir->getPar(0, 0, 3);
0121     std::cout << "propagate to lay=" << ilay << " start from x=" << fir->getPar(0, 0, 0)
0122               << " y=" << fir->getPar(0, 0, 1) << " z=" << fir->getPar(0, 0, 2)
0123               << " r=" << getHypot(fir->getPar(0, 0, 0), fir->getPar(0, 0, 1))
0124               << " px=" << pt * std::cos(fir->getPar(0, 0, 4)) << " py=" << pt * std::sin(fir->getPar(0, 0, 4))
0125               << " pz=" << pt / std::tan(fir->getPar(0, 0, 5)) << " pT=" << pt << std::endl;
0126   }
0127
0128   void post_prop_print(int ilay, MkBase *fir) {
0129     std::cout << "propagate to lay=" << ilay << " arrive at x=" << fir->getPar(0, 1, 0) << " y=" << fir->getPar(0, 1, 1)
0130               << " z=" << fir->getPar(0, 1, 2) << " r=" << getHypot(fir->getPar(0, 1, 0), fir->getPar(0, 1, 1))
0131               << std::endl;
0132   }
0133
0134   void print_seed(const Track &seed) {
0135     std::cout << "MX - found seed with label=" << seed.label() << " nHits=" << seed.nFoundHits()
0136               << " chi2=" << seed.chi2() << " posEta=" << seed.posEta() << " posPhi=" << seed.posPhi()
0137               << " posR=" << seed.posR() << " posZ=" << seed.z() << " pT=" << seed.pT() << std::endl;
0138   }
0139
0140   void print_seed2(const TrackCand &seed) {
0141     std::cout << "MX - found seed with nFoundHits=" << seed.nFoundHits() << " chi2=" << seed.chi2() << " x=" << seed.x()
0142               << " y=" << seed.y() << " z=" << seed.z() << " px=" << seed.px() << " py=" << seed.py()
0143               << " pz=" << seed.pz() << " pT=" << seed.pT() << std::endl;
0144   }
0145
0146   void print_seeds(const TrackVec &seeds) {
0147     std::cout << "found total seeds=" << seeds.size() << std::endl;
0148     for (auto &&seed : seeds) {
0149       print_seed(seed);
0150     }
0151   }
0152
0153   [[maybe_unused]] void print_seeds(const EventOfCombCandidates &event_of_comb_cands) {
0154     for (int iseed = 0; iseed < event_of_comb_cands.size(); iseed++) {
0155       print_seed2(event_of_comb_cands[iseed].front());
0156     }
0157   }
0158 #endif
0159
0160   bool sortCandByScore(const TrackCand &cand1, const TrackCand &cand2) {
0161     return mkfit::sortByScoreTrackCand(cand1, cand2);
0162   }
0163
0164 #ifdef RNT_DUMP_MkF_SelHitIdcs
0165   constexpr bool alwaysUseHitSelectionV2 = true;
0166 #else
0167   constexpr bool alwaysUseHitSelectionV2 = false;
0168 #endif
0169 }  // end unnamed namespace
0170
0171 //------------------------------------------------------------------------------
0172 // Constructor and destructor
0173 //------------------------------------------------------------------------------
0174
0175 namespace mkfit {
0176
0177   std::unique_ptr<MkBuilder> MkBuilder::make_builder(bool silent) { return std::make_unique<MkBuilder>(silent); }
0178
0179   void MkBuilder::populate() { g_exe_ctx.populate(Config::numThreadsFinder); }
0180   void MkBuilder::clear() { g_exe_ctx.clear(); }
0181
0182   std::pair<int, int> MkBuilder::max_hits_layer(const EventOfHits &eoh) const {
0183     int maxN = 0;
0184     int maxL = 0;
0185     for (int l = 0; l < eoh.nLayers(); ++l) {
0186       int lsize = eoh[l].nHits();
0187       if (lsize > maxN) {
0188         maxN = lsize;
0189         maxL = eoh[l].layer_id();
0190       }
0191     }
0192     return {maxN, maxL};
0193   }
0194
0195   int MkBuilder::total_cands() const {
0196     int res = 0;
0197     for (int i = 0; i < m_event_of_comb_cands.size(); ++i)
0198       res += m_event_of_comb_cands[i].size();
0199     return res;
0200   }
0201
0202   //------------------------------------------------------------------------------
0203   // Common functions
0204   //------------------------------------------------------------------------------
0205
0206   void MkBuilder::begin_event(MkJob *job, Event *ev, const char *build_type) {
0207     m_nan_n_silly_per_layer_count = 0;
0208
0209     m_job = job;
0210     m_event = ev;
0211
0212     m_seedEtaSeparators.resize(m_job->num_regions());
0213     m_seedMinLastLayer.resize(m_job->num_regions());
0214     m_seedMaxLastLayer.resize(m_job->num_regions());
0215
0216     for (int i = 0; i < m_job->num_regions(); ++i) {
0217       m_seedEtaSeparators[i] = 0;
0218       m_seedMinLastLayer[i] = 9999;
0219       m_seedMaxLastLayer[i] = 0;
0220     }
0221
0222     if (!m_silent) {
0223       std::cout << "MkBuilder building tracks with '" << build_type << "'"
0224                 << ", iteration_index=" << job->m_iter_config.m_iteration_index
0225                 << ", track_algorithm=" << job->m_iter_config.m_track_algorithm << std::endl;
0226     }
0227   }
0228
0229   void MkBuilder::end_event() {
0230     m_job = nullptr;
0231     m_event = nullptr;
0232   }
0233
0234   void MkBuilder::release_memory() {
0235     TrackVec tmp;
0236     m_tracks.swap(tmp);
0237     m_event_of_comb_cands.releaseMemory();
0238   }
0239
0240   void MkBuilder::import_seeds(const TrackVec &in_seeds,
0241                                const bool seeds_sorted,
0242                                std::function<insert_seed_foo> insert_seed) {
0243     // bool debug = true;
0244
0245     const int size = in_seeds.size();
0246
0247     IterationSeedPartition part(size);
0248     std::vector<unsigned> ranks;
0249     if (!seeds_sorted) {
0250       // We don't care about bins in phi, use low N to reduce overall number of bins.
0251       axis_pow2_u1<float, unsigned short, 10, 4> ax_phi(-Const::PI, Const::PI);
0252       axis<float, unsigned short, 8, 8> ax_eta(-3.0, 3.0, 64u);
0253       binnor<unsigned int, decltype(ax_phi), decltype(ax_eta), 20, 12> phi_eta_binnor(ax_phi, ax_eta);
0254       part.m_phi_eta_foo = [&](float phi, float eta) { phi_eta_binnor.register_entry_safe(phi, eta); };
0255
0256       phi_eta_binnor.begin_registration(size);
0257       m_job->m_iter_config.m_seed_partitioner(m_job->m_trk_info, in_seeds, m_job->m_event_of_hits, part);
0258       phi_eta_binnor.finalize_registration();
0259       ranks.swap(phi_eta_binnor.m_ranks);
0260     } else {
0261       m_job->m_iter_config.m_seed_partitioner(m_job->m_trk_info, in_seeds, m_job->m_event_of_hits, part);
0262     }
0263
0264     for (int i = 0; i < size; ++i) {
0265       int j = seeds_sorted ? i : ranks[i];
0266       int reg = part.m_region[j];
0267       ++m_seedEtaSeparators[reg];
0268     }
0269
0270     // Sum up region counts to contain actual ending indices and prepare insertion cursors.
0271     // Fix min/max layers.
0272     std::vector<int> seed_cursors(m_job->num_regions());
0273     for (int reg = 1; reg < m_job->num_regions(); ++reg) {
0274       seed_cursors[reg] = m_seedEtaSeparators[reg - 1];
0275       m_seedEtaSeparators[reg] += m_seedEtaSeparators[reg - 1];
0276     }
0277
0278     // Actually imports seeds, detect last-hit layer range per region.
0279     for (int i = 0; i < size; ++i) {
0280       int j = seeds_sorted ? i : ranks[i];
0281       int reg = part.m_region[j];
0282       const Track &seed = in_seeds[j];
0283       insert_seed(seed, j, reg, seed_cursors[reg]++);
0284
0285       HitOnTrack hot = seed.getLastHitOnTrack();
0286       m_seedMinLastLayer[reg] = std::min(m_seedMinLastLayer[reg], hot.layer);
0287       m_seedMaxLastLayer[reg] = std::max(m_seedMaxLastLayer[reg], hot.layer);
0288     }
0289
0290     // Fix min/max layers
0291     for (int reg = 0; reg < m_job->num_regions(); ++reg) {
0292       if (m_seedMinLastLayer[reg] == 9999)
0293         m_seedMinLastLayer[reg] = -1;
0294       if (m_seedMaxLastLayer[reg] == 0)
0295         m_seedMaxLastLayer[reg] = -1;
0296     }
0297
0298     // clang-format off
0299     dprintf("MkBuilder::import_seeds finished import of %d seeds (last seeding layer min, max):\n"
0300             "  ec- = %d(%d,%d), t- = %d(%d,%d), brl = %d(%d,%d), t+ = %d(%d,%d), ec+ = %d(%d,%d).\n",
0301             size,
0302             m_seedEtaSeparators[0],                          m_seedMinLastLayer[0], m_seedMaxLastLayer[0],
0303             m_seedEtaSeparators[1] - m_seedEtaSeparators[0], m_seedMinLastLayer[1], m_seedMaxLastLayer[1],
0304             m_seedEtaSeparators[2] - m_seedEtaSeparators[1], m_seedMinLastLayer[2], m_seedMaxLastLayer[2],
0305             m_seedEtaSeparators[3] - m_seedEtaSeparators[2], m_seedMinLastLayer[3], m_seedMaxLastLayer[3],
0306             m_seedEtaSeparators[4] - m_seedEtaSeparators[3], m_seedMinLastLayer[4], m_seedMaxLastLayer[4]);
0307     // dcall(print_seeds(m_event_of_comb_cands));
0308     // clang-format on
0309   }
0310
0311   //------------------------------------------------------------------------------
0312
0313   int MkBuilder::filter_comb_cands(filter_candidates_func filter, bool attempt_all_cands) {
0314     EventOfCombCandidates &eoccs = m_event_of_comb_cands;
0315     int i = 0, place_pos = 0;
0316
0317     dprintf("MkBuilder::filter_comb_cands Entering filter size eoccs.size=%d\n", eoccs.size());
0318
0319     std::vector<int> removed_cnts(m_job->num_regions());
0320     while (i < eoccs.size()) {
0321       if (eoccs.cands_in_backward_rep())
0322         eoccs[i].repackCandPostBkwSearch(0);
0323       bool passed = filter(eoccs[i].front(), *m_job);
0324
0325       if (!passed && attempt_all_cands) {
0326         for (int j = 1; j < (int)eoccs[i].size(); ++j) {
0327           if (eoccs.cands_in_backward_rep())
0328             eoccs[i].repackCandPostBkwSearch(j);
0329           if (filter(eoccs[i][j], *m_job)) {
0330             eoccs[i][0] = eoccs[i][j];  // overwrite front, no need to std::swap() them
0331             passed = true;
0332             break;
0333           }
0334         }
0335       }
0336       if (passed) {
0337         if (place_pos != i)
0338           std::swap(eoccs[place_pos], eoccs[i]);
0339         ++place_pos;
0340       } else {
0341         assert(eoccs[i].front().getEtaRegion() < m_job->num_regions());
0342         ++removed_cnts[eoccs[i].front().getEtaRegion()];
0343       }
0344       ++i;
0345     }
0346
0347     int n_removed = 0;
0348     for (int reg = 0; reg < m_job->num_regions(); ++reg) {
0349       dprintf("MkBuilder::filter_comb_cands reg=%d: n_rem_was=%d removed_in_r=%d n_rem=%d, es_was=%d es_new=%d\n",
0350               reg,
0351               n_removed,
0352               removed_cnts[reg],
0353               n_removed + removed_cnts[reg],
0354               m_seedEtaSeparators[reg],
0355               m_seedEtaSeparators[reg] - n_removed - removed_cnts[reg]);
0356
0357       n_removed += removed_cnts[reg];
0358       m_seedEtaSeparators[reg] -= n_removed;
0359     }
0360
0361     eoccs.resizeAfterFiltering(n_removed);
0362
0363     dprintf("MkBuilder::filter_comb_cands n_removed = %d, eoccs.size=%d\n", n_removed, eoccs.size());
0364
0365     return n_removed;
0366   }
0367
0368   void MkBuilder::find_min_max_hots_size() {
0369     const EventOfCombCandidates &eoccs = m_event_of_comb_cands;
0370     int min[5], max[5], gmin = 0, gmax = 0;
0371     int i = 0;
0372     for (int reg = 0; reg < 5; ++reg) {
0373       min[reg] = 9999;
0374       max[reg] = 0;
0375       for (; i < m_seedEtaSeparators[reg]; i++) {
0376         min[reg] = std::min(min[reg], eoccs[i].hotsSize());
0377         max[reg] = std::max(max[reg], eoccs[i].hotsSize());
0378       }
0379       gmin = std::max(gmin, min[reg]);
0380       gmax = std::max(gmax, max[reg]);
0381     }
0382     // clang-format off
0383     printf("MkBuilder::find_min_max_hots_size MIN %3d -- [ %3d | %3d | %3d | %3d | %3d ]   "
0384            "MAX %3d -- [ %3d | %3d | %3d | %3d | %3d ]\n",
0385            gmin, min[0], min[1], min[2], min[3], min[4],
0386            gmax, max[0], max[1], max[2], max[3], max[4]);
0387     // clang-format on
0388   }
0389
0390   void MkBuilder::select_best_comb_cands(bool clear_m_tracks, bool remove_missing_hits) {
0391     if (clear_m_tracks)
0392       m_tracks.clear();
0393     export_best_comb_cands(m_tracks, remove_missing_hits);
0394   }
0395
0396   void MkBuilder::export_best_comb_cands(TrackVec &out_vec, bool remove_missing_hits) {
0397     const EventOfCombCandidates &eoccs = m_event_of_comb_cands;
0398     out_vec.reserve(out_vec.size() + eoccs.size());
0399     for (int i = 0; i < eoccs.size(); i++) {
0400       // Take the first candidate, if it exists.
0401       if (!eoccs[i].empty()) {
0402         const TrackCand &bcand = eoccs[i].front();
0403         out_vec.emplace_back(bcand.exportTrack(remove_missing_hits));
0404       }
0405     }
0406   }
0407
0408   void MkBuilder::export_tracks(TrackVec &out_vec) {
0409     out_vec.reserve(out_vec.size() + m_tracks.size());
0410     for (auto &t : m_tracks) {
0411       out_vec.emplace_back(t);
0412     }
0413   }
0414
0415   //------------------------------------------------------------------------------
0416   // PrepareSeeds
0417   //------------------------------------------------------------------------------
0418
0419   void MkBuilder::seed_post_cleaning(TrackVec &tv) {
0420     if constexpr (Const::nan_n_silly_check_seeds) {
0421       int count = 0;
0422
0423       for (int i = 0; i < (int)tv.size(); ++i) {
0424         bool silly = tv[i].hasSillyValues(Const::nan_n_silly_print_bad_seeds,
0425                                           Const::nan_n_silly_fixup_bad_seeds,
0426                                           "Post-cleaning seed silly value check and fix");
0427         if (silly) {
0428           ++count;
0429           if constexpr (Const::nan_n_silly_remove_bad_seeds) {
0430             // XXXX MT
0431             // Could do somethin smarter here: set as Stopped ?  check in seed cleaning ?
0432             tv.erase(tv.begin() + i);
0433             --i;
0434           }
0435         }
0436       }
0437
0438       if (count > 0 && !m_silent) {
0439         printf("Nan'n'Silly detected %d silly seeds (fix=%d, remove=%d).\n",
0440                count,
0441                Const::nan_n_silly_fixup_bad_seeds,
0442                Const::nan_n_silly_remove_bad_seeds);
0443       }
0444     }
0445   }
0446
0447   //------------------------------------------------------------------------------
0448   // FindTracksBestHit
0449   //------------------------------------------------------------------------------
0450
0451   void MkBuilder::find_tracks_load_seeds_BH(const TrackVec &in_seeds, const bool seeds_sorted) {
0452     // bool debug = true;
0453     assert(!in_seeds.empty());
0454     m_tracks.resize(in_seeds.size());
0455
0456     import_seeds(in_seeds, seeds_sorted, [&](const Track &seed, int seed_idx, int region, int pos) {
0457       m_tracks[pos] = seed;
0458       m_tracks[pos].setNSeedHits(seed.nTotalHits());
0459       m_tracks[pos].setEtaRegion(region);
0460     });
0461
0462     //dump seeds
0463     dcall(print_seeds(m_tracks));
0464   }
0465
0466   void MkBuilder::findTracksBestHit(SteeringParams::IterationType_e iteration_dir) {
0467     // bool debug = true;
0468
0469     TrackVec &cands = m_tracks;
0470
0471     tbb::parallel_for_each(m_job->regions_begin(), m_job->regions_end(), [&](int region) {
0472       if (iteration_dir == SteeringParams::IT_BkwSearch && !m_job->steering_params(region).has_bksearch_plan()) {
0473         printf("No backward search plan for region %d\n", region);
0474         return;
0475       }
0476
0477       // XXXXXX Select endcap / barrel only ...
0478       // if (region != TrackerInfo::Reg_Endcap_Neg && region != TrackerInfo::Reg_Endcap_Pos)
0479       // if (region != TrackerInfo::Reg_Barrel)
0480       //   return;
0481
0482       const SteeringParams &st_par = m_job->steering_params(region);
0483       const TrackerInfo &trk_info = m_job->m_trk_info;
0484       const PropagationConfig &prop_config = trk_info.prop_config();
0485
0486       const RegionOfSeedIndices rosi(m_seedEtaSeparators, region);
0487
0488       tbb::parallel_for(rosi.tbb_blk_rng_vec(), [&](const tbb::blocked_range<int> &blk_rng) {
0489         auto mkfndr = g_exe_ctx.m_finders.makeOrGet();
0490
0491         RangeOfSeedIndices rng = rosi.seed_rng(blk_rng);
0492
0493         std::vector<int> trk_idcs(NN);  // track indices in Matriplex
0494         std::vector<int> trk_llay(NN);  // last layer on input track
0495
0496         while (rng.valid()) {
0497           dprint(std::endl << "processing track=" << rng.m_beg << ", label=" << cands[rng.m_beg].label());
0498
0499           int prev_layer = 9999;
0500
0501           for (int i = rng.m_beg, ii = 0; i < rng.m_end; ++i, ++ii) {
0502             int llay = cands[i].getLastHitLyr();
0503             trk_llay[ii] = llay;
0504             prev_layer = std::min(prev_layer, llay);
0505
0506             dprintf("  %2d %2d %2d lay=%3d prev_layer=%d\n", ii, i, cands[i].label(), llay, prev_layer);
0507           }
0508           int curr_tridx = 0;
0509
0510           auto layer_plan_it = st_par.make_iterator(iteration_dir);
0511
0512           dprintf("Made iterator for %d, first layer=%d ... end layer=%d\n",
0513                   iteration_dir,
0514                   layer_plan_it.layer(),
0515                   layer_plan_it.last_layer());
0516
0517           assert(layer_plan_it.is_pickup_only());
0518
0519           int curr_layer = layer_plan_it.layer();
0520
0521           mkfndr->m_Stopped.setVal(0);
0522
0523           // Loop over layers, starting from after the seed.
0524           // Consider inverting loop order and make layer outer, need to
0525           // trade off hit prefetching with copy-out of candidates.
0526           while (++layer_plan_it) {
0527             prev_layer = curr_layer;
0528             curr_layer = layer_plan_it.layer();
0529             mkfndr->setup(prop_config,
0530                           m_job->m_iter_config,
0531                           m_job->m_iter_config.m_params,
0532                           m_job->m_iter_config.m_layer_configs[curr_layer],
0533                           st_par,
0534                           m_job->get_mask_for_layer(curr_layer),
0535                           m_event,
0536                           region,
0537                           m_job->m_in_fwd);
0538
0539             const LayerOfHits &layer_of_hits = m_job->m_event_of_hits[curr_layer];
0540             const LayerInfo &layer_info = trk_info.layer(curr_layer);
0541             const FindingFoos &fnd_foos = FindingFoos::get_finding_foos(layer_info.is_barrel());
0542             dprint("at layer " << curr_layer << ", nHits in layer " << layer_of_hits.nHits());
0543
0544             // Pick up seeds that become active on current layer -- unless already fully loaded.
0545             if (curr_tridx < rng.n_proc()) {
0546               int prev_tridx = curr_tridx;
0547
0548               for (int i = rng.m_beg, ii = 0; i < rng.m_end; ++i, ++ii) {
0549                 if (trk_llay[ii] == prev_layer)
0550                   trk_idcs[curr_tridx++] = i;
0551               }
0552               if (curr_tridx > prev_tridx) {
0553                 dprintf("added %d seeds, started with %d\n", curr_tridx - prev_tridx, prev_tridx);
0554
0555                 mkfndr->inputTracksAndHitIdx(cands, trk_idcs, prev_tridx, curr_tridx, false, prev_tridx);
0556               }
0557             }
0558
0559             if (layer_plan_it.is_pickup_only())
0560               continue;
0561
0562             dcall(pre_prop_print(curr_layer, mkfndr.get()));
0563
0564             mkfndr->clearFailFlag();
0565             (mkfndr.get()->*fnd_foos.m_propagate_foo)(
0566                 layer_info.propagate_to(), curr_tridx, prop_config.finding_inter_layer_pflags);
0567
0568             dcall(post_prop_print(curr_layer, mkfndr.get()));
0569
0570             mkfndr->selectHitIndices(layer_of_hits, curr_tridx);
0571
0572             // Stop low-pT tracks that can not reach the current barrel layer.
0573             if (layer_info.is_barrel()) {
0574               const float r_min_sqr = layer_info.rin() * layer_info.rin();
0575               for (int i = 0; i < curr_tridx; ++i) {
0576                 if (!mkfndr->m_Stopped[i]) {
0577                   if (mkfndr->radiusSqr(i, MkBase::iP) < r_min_sqr) {
0578                     if (region == TrackerInfo::Reg_Barrel) {
0579                       mkfndr->m_Stopped[i] = 1;
0580                       mkfndr->outputTrackAndHitIdx(cands[rng.m_beg + i], i, false);
0581                     }
0582                     mkfndr->m_XWsrResult[i].m_wsr = WSR_Outside;
0583                     mkfndr->m_XHitSize[i] = 0;
0584                   }
0585                 } else {  // make sure we don't add extra work for AddBestHit
0586                   mkfndr->m_XWsrResult[i].m_wsr = WSR_Outside;
0587                   mkfndr->m_XHitSize[i] = 0;
0588                 }
0589               }
0590             }
0591
0592             // make candidates with best hit
0593             dprint("make new candidates");
0594
0595             mkfndr->addBestHit(layer_of_hits, curr_tridx, fnd_foos);
0596
0597             // Stop tracks that have reached N_max_holes.
0598             for (int i = 0; i < curr_tridx; ++i) {
0599               if (!mkfndr->m_Stopped[i] && mkfndr->bestHitLastHoT(i).index == -2) {
0600                 mkfndr->m_Stopped[i] = 1;
0601                 mkfndr->outputTrackAndHitIdx(cands[rng.m_beg + i], i, false);
0602               }
0603             }
0604
0605           }  // end of layer loop
0606
0607           mkfndr->outputNonStoppedTracksAndHitIdx(cands, trk_idcs, 0, curr_tridx, false);
0608
0609           ++rng;
0610         }  // end of loop over candidates in a tbb chunk
0611
0612         mkfndr->release();
0613       });  // end parallel_for over candidates in a region
0614     });    // end of parallel_for_each over regions
0615   }
0616
0617   //------------------------------------------------------------------------------
0618   // FindTracksCombinatorial: Standard TBB and CloneEngine TBB
0619   //------------------------------------------------------------------------------
0620
0621   void MkBuilder::find_tracks_load_seeds(const TrackVec &in_seeds, const bool seeds_sorted) {
0622     // This will sort seeds according to iteration configuration.
0623     assert(!in_seeds.empty());
0624     m_tracks.clear();  // m_tracks can be used for BkFit.
0625
0626     m_event_of_comb_cands.reset((int)in_seeds.size(), m_job->max_max_cands());
0627
0628     import_seeds(in_seeds, seeds_sorted, [&](const Track &seed, int seed_idx, int region, int pos) {
0629       m_event_of_comb_cands.insertSeed(seed, seed_idx, m_job->steering_params(region).m_track_scorer, region, pos);
0630     });
0631   }
0632
0633   int MkBuilder::find_tracks_unroll_candidates(std::vector<std::pair<int, int>> &seed_cand_vec,
0634                                                int start_seed,
0635                                                int end_seed,
0636                                                int layer,
0637                                                int prev_layer,
0638                                                bool pickup_only,
0639                                                SteeringParams::IterationType_e iteration_dir) {
0640     int silly_count = 0;
0641
0642     seed_cand_vec.clear();
0643
0644     auto &iter_params = (iteration_dir == SteeringParams::IT_BkwSearch) ? m_job->m_iter_config.m_backward_params
0645                                                                         : m_job->m_iter_config.m_params;
0646
0647     for (int iseed = start_seed; iseed < end_seed; ++iseed) {
0648       CombCandidate &ccand = m_event_of_comb_cands[iseed];
0649
0650       if (ccand.state() == CombCandidate::Dormant && ccand.pickupLayer() == prev_layer) {
0651         ccand.setState(CombCandidate::Finding);
0652       }
0653       if (!pickup_only && ccand.state() == CombCandidate::Finding) {
0654         bool active = false;
0655         for (int ic = 0; ic < (int)ccand.size(); ++ic) {
0656           if (ccand[ic].getLastHitIdx() != -2) {
0657             // Stop candidates with pT<X GeV
0658             if (ccand[ic].pT() < iter_params.minPtCut) {
0659               ccand[ic].addHitIdx(-2, layer, 0.0f);
0660               continue;
0661             }
0662             // Check if the candidate is close to it's max_r, pi/2 - 0.2 rad (11.5 deg)
0663             if (iteration_dir == SteeringParams::IT_FwdSearch && ccand[ic].pT() < 1.2) {
0664               const float dphi = std::abs(ccand[ic].posPhi() - ccand[ic].momPhi());
0665               if (ccand[ic].posRsq() > 625.f && dphi > 1.371f && dphi < 4.512f) {
0666                 // printf("Stopping cand at r=%f, posPhi=%.1f momPhi=%.2f pt=%.2f emomEta=%.2f\n",
0667                 //        ccand[ic].posR(), ccand[ic].posPhi(), ccand[ic].momPhi(), ccand[ic].pT(), ccand[ic].momEta());
0668                 ccand[ic].addHitIdx(-2, layer, 0.0f);
0669                 continue;
0670               }
0671             }
0672
0673             active = true;
0674             seed_cand_vec.push_back(std::pair<int, int>(iseed, ic));
0675             ccand[ic].resetOverlaps();
0676
0677             if constexpr (Const::nan_n_silly_check_cands_every_layer) {
0678               if (ccand[ic].hasSillyValues(Const::nan_n_silly_print_bad_cands_every_layer,
0679                                            Const::nan_n_silly_fixup_bad_cands_every_layer,
0680                                            "Per layer silly check"))
0681                 ++silly_count;
0682             }
0683           }
0684         }
0685         if (!active) {
0686           ccand.setState(CombCandidate::Finished);
0687         }
0688       }
0689     }
0690
0691     if constexpr (Const::nan_n_silly_check_cands_every_layer && silly_count > 0) {
0692       m_nan_n_silly_per_layer_count += silly_count;
0693     }
0694
0695     return seed_cand_vec.size();
0696   }
0697
0698   void MkBuilder::find_tracks_handle_missed_layers(MkFinder *mkfndr,
0699                                                    const LayerInfo &layer_info,
0700                                                    std::vector<std::vector<TrackCand>> &tmp_cands,
0701                                                    const std::vector<std::pair<int, int>> &seed_cand_idx,
0702                                                    const int region,
0703                                                    const int start_seed,
0704                                                    const int itrack,
0705                                                    const int end) {
0706     // XXXX-1 If I miss a layer, insert the original track into tmp_cands
0707     // AND do not do it in FindCandidates as the position can be badly
0708     // screwed by then. See comment there, too.
0709     // One could also do a pre-check ... so as not to use up a slot.
0710
0711     // bool debug = true;
0712
0713     for (int ti = itrack; ti < end; ++ti) {
0714       TrackCand &cand = m_event_of_comb_cands[seed_cand_idx[ti].first][seed_cand_idx[ti].second];
0715       WSR_Result &w = mkfndr->m_XWsrResult[ti - itrack];
0716
0717       // Low pT tracks can miss a barrel layer ... and should be stopped
0718       dprintf("WSR Check label %d, seed %d, cand %d score %f -> wsr %d, in_gap %d\n",
0719               cand.label(),
0720               seed_cand_idx[ti].first,
0721               seed_cand_idx[ti].second,
0722               cand.score(),
0723               w.m_wsr,
0724               w.m_in_gap);
0725
0726       if (w.m_wsr == WSR_Failed) {
0727         // Fake outside so it does not get processed in FindTracks BH/Std/CE.
0728         // [ Should add handling of WSR_Failed there, perhaps. ]
0729         w.m_wsr = WSR_Outside;
0730
0731         if (layer_info.is_barrel()) {
0732           dprintf("Barrel cand propagation failed, got to r=%f ... layer is %f - %f\n",
0733                   mkfndr->radius(ti - itrack, MkBase::iP),
0734                   layer_info.rin(),
0735                   layer_info.rout());
0736           // In barrel region, create a stopped replica. In transition region keep the original copy
0737           // as there is still a chance to hit endcaps.
0738           tmp_cands[seed_cand_idx[ti].first - start_seed].push_back(cand);
0739           if (region == TrackerInfo::Reg_Barrel) {
0740             dprintf(" creating extra stopped held back candidate\n");
0741             tmp_cands[seed_cand_idx[ti].first - start_seed].back().addHitIdx(-2, layer_info.layer_id(), 0);
0742           }
0743         }
0744         // Never happens for endcap / propToZ
0745       } else if (w.m_wsr == WSR_Outside) {
0746         dprintf(" creating extra held back candidate\n");
0747         tmp_cands[seed_cand_idx[ti].first - start_seed].push_back(cand);
0748       } else if (w.m_wsr == WSR_Edge) {
0749         // Do nothing special here, this case is handled also in MkFinder:findTracks()
0750       }
0751     }
0752   }
0753
0754   //------------------------------------------------------------------------------
0755   // FindTracksCombinatorial: Standard TBB
0756   //------------------------------------------------------------------------------
0757
0758   void MkBuilder::findTracksStandard(SteeringParams::IterationType_e iteration_dir) {
0759     // debug = true;
0760
0761     EventOfCombCandidates &eoccs = m_event_of_comb_cands;
0762
0763     tbb::parallel_for_each(m_job->regions_begin(), m_job->regions_end(), [&](int region) {
0764       if (iteration_dir == SteeringParams::IT_BkwSearch && !m_job->steering_params(region).has_bksearch_plan()) {
0765         printf("No backward search plan for region %d\n", region);
0766         return;
0767       }
0768
0769       const TrackerInfo &trk_info = m_job->m_trk_info;
0770       const SteeringParams &st_par = m_job->steering_params(region);
0771       const IterationParams &params = m_job->params();
0772       const PropagationConfig &prop_config = trk_info.prop_config();
0773
0774       const RegionOfSeedIndices rosi(m_seedEtaSeparators, region);
0775
0776       // adaptive seeds per task based on the total estimated amount of work to divide among all threads
0777       const int adaptiveSPT = std::clamp(
0778           Config::numThreadsEvents * eoccs.size() / Config::numThreadsFinder + 1, 4, Config::numSeedsPerTask);
0779       dprint("adaptiveSPT " << adaptiveSPT << " fill " << rosi.count() << "/" << eoccs.size() << " region " << region);
0780
0781       // loop over seeds
0782       tbb::parallel_for(rosi.tbb_blk_rng_std(adaptiveSPT), [&](const tbb::blocked_range<int> &seeds) {
0783         auto mkfndr = g_exe_ctx.m_finders.makeOrGet();
0784
0785         const int start_seed = seeds.begin();
0786         const int end_seed = seeds.end();
0787         const int n_seeds = end_seed - start_seed;
0788
0789         std::vector<std::vector<TrackCand>> tmp_cands(n_seeds);
0790         for (size_t iseed = 0; iseed < tmp_cands.size(); ++iseed) {
0791           tmp_cands[iseed].reserve(2 * params.maxCandsPerSeed);  //factor 2 seems reasonable to start with
0792         }
0793
0794         std::vector<std::pair<int, int>> seed_cand_idx;
0795         seed_cand_idx.reserve(n_seeds * params.maxCandsPerSeed);
0796
0797         auto layer_plan_it = st_par.make_iterator(iteration_dir);
0798
0799         dprintf("Made iterator for %d, first layer=%d ... end layer=%d\n",
0800                 iteration_dir,
0801                 layer_plan_it.layer(),
0802                 layer_plan_it.last_layer());
0803
0804         assert(layer_plan_it.is_pickup_only());
0805
0806         int curr_layer = layer_plan_it.layer(), prev_layer;
0807
0808         dprintf("\nMkBuilder::FindTracksStandard region=%d, seed_pickup_layer=%d, first_layer=%d\n",
0809                 region,
0810                 curr_layer,
0811                 layer_plan_it.next_layer());
0812
0813         auto &iter_params = (iteration_dir == SteeringParams::IT_BkwSearch) ? m_job->m_iter_config.m_backward_params
0814                                                                             : m_job->m_iter_config.m_params;
0815
0816         // Loop over layers, starting from after the seed.
0817         while (++layer_plan_it) {
0818           prev_layer = curr_layer;
0819           curr_layer = layer_plan_it.layer();
0820           mkfndr->setup(prop_config,
0821                         m_job->m_iter_config,
0822                         iter_params,
0823                         m_job->m_iter_config.m_layer_configs[curr_layer],
0824                         st_par,
0825                         m_job->get_mask_for_layer(curr_layer),
0826                         m_event,
0827                         region,
0828                         m_job->m_in_fwd);
0829
0830           const LayerOfHits &layer_of_hits = m_job->m_event_of_hits[curr_layer];
0831           const LayerInfo &layer_info = trk_info.layer(curr_layer);
0832           const FindingFoos &fnd_foos = FindingFoos::get_finding_foos(layer_info.is_barrel());
0833
0834           dprintf("\n* Processing layer %d\n", curr_layer);
0835           mkfndr->begin_layer(layer_of_hits);
0836
0837           int theEndCand = find_tracks_unroll_candidates(seed_cand_idx,
0838                                                          start_seed,
0839                                                          end_seed,
0840                                                          curr_layer,
0841                                                          prev_layer,
0842                                                          layer_plan_it.is_pickup_only(),
0843                                                          iteration_dir);
0844
0845           dprintf("  Number of candidates to process: %d, nHits in layer: %d\n", theEndCand, layer_of_hits.nHits());
0846
0847           if (layer_plan_it.is_pickup_only() || theEndCand == 0)
0848             continue;
0849
0850           // vectorized loop
0851           for (int itrack = 0; itrack < theEndCand; itrack += NN) {
0852             int end = std::min(itrack + NN, theEndCand);
0853
0854             dprint("processing track=" << itrack << ", label="
0855                                        << eoccs[seed_cand_idx[itrack].first][seed_cand_idx[itrack].second].label());
0856
0857             //fixme find a way to deal only with the candidates needed in this thread
0858             mkfndr->inputTracksAndHitIdx(eoccs.refCandidates(), seed_cand_idx, itrack, end, false);
0859
0860             //propagate to layer
0861             dcall(pre_prop_print(curr_layer, mkfndr.get()));
0862
0863             mkfndr->clearFailFlag();
0864             (mkfndr.get()->*fnd_foos.m_propagate_foo)(
0865                 layer_info.propagate_to(), end - itrack, prop_config.finding_inter_layer_pflags);
0866
0867             dcall(post_prop_print(curr_layer, mkfndr.get()));
0868
0869             dprint("now get hit range");
0870
0871             if (alwaysUseHitSelectionV2 || iter_params.useHitSelectionV2)
0872               mkfndr->selectHitIndicesV2(layer_of_hits, end - itrack);
0873             else
0874               mkfndr->selectHitIndices(layer_of_hits, end - itrack);
0875
0876             find_tracks_handle_missed_layers(
0877                 mkfndr.get(), layer_info, tmp_cands, seed_cand_idx, region, start_seed, itrack, end);
0878
0879             dprint("make new candidates");
0880             mkfndr->findCandidates(layer_of_hits, tmp_cands, start_seed, end - itrack, fnd_foos);
0881
0882           }  //end of vectorized loop
0883
0884           // sort the input candidates
0885           for (int is = 0; is < n_seeds; ++is) {
0886             dprint("dump seed n " << is << " with N_input_candidates=" << tmp_cands[is].size());
0887
0888             std::sort(tmp_cands[is].begin(), tmp_cands[is].end(), sortCandByScore);
0889           }
0890
0891           // now fill out the output candidates
0892           for (int is = 0; is < n_seeds; ++is) {
0893             if (!tmp_cands[is].empty()) {
0894               eoccs[start_seed + is].clear();
0895
0896               // Put good candidates into eoccs, process -2 candidates.
0897               int n_placed = 0;
0898               bool first_short = true;
0899               for (int ii = 0; ii < (int)tmp_cands[is].size() && n_placed < params.maxCandsPerSeed; ++ii) {
0900                 TrackCand &tc = tmp_cands[is][ii];
0901
0902                 // See if we have an overlap hit available, but only if we have a true hit in this layer
0903                 // and pT is above the pTCutOverlap
0904                 if (tc.pT() > params.pTCutOverlap && tc.getLastHitLyr() == curr_layer && tc.getLastHitIdx() >= 0) {
0905                   CombCandidate &ccand = eoccs[start_seed + is];
0906
0907                   HitMatch *hm = ccand[tc.originIndex()].findOverlap(
0908                       tc.getLastHitIdx(), layer_of_hits.refHit(tc.getLastHitIdx()).detIDinLayer());
0909
0910                   if (hm) {
0911                     tc.addHitIdx(hm->m_hit_idx, curr_layer, hm->m_chi2);
0912                     tc.incOverlapCount();
0913                   }
0914                 }
0915
0916                 if (tc.getLastHitIdx() != -2) {
0917                   eoccs[start_seed + is].emplace_back(tc);
0918                   ++n_placed;
0919                 } else if (first_short) {
0920                   first_short = false;
0921                   if (tc.score() > eoccs[start_seed + is].refBestShortCand().score()) {
0922                     eoccs[start_seed + is].setBestShortCand(tc);
0923                   }
0924                 }
0925               }
0926
0927               tmp_cands[is].clear();
0928             }
0929           }
0930           mkfndr->end_layer();
0931         }  // end of layer loop
0932         mkfndr->release();
0933
0934         // final sorting
0935         for (int iseed = start_seed; iseed < end_seed; ++iseed) {
0936           eoccs[iseed].mergeCandsAndBestShortOne(m_job->params(), st_par.m_track_scorer, true, true);
0937         }
0938       });  // end parallel-for over chunk of seeds within region
0939     });    // end of parallel-for-each over eta regions
0940
0941     // debug = false;
0942   }
0943
0944   //------------------------------------------------------------------------------
0945   // FindTracksCombinatorial: CloneEngine TBB
0946   //------------------------------------------------------------------------------
0947
0948   void MkBuilder::findTracksCloneEngine(SteeringParams::IterationType_e iteration_dir) {
0949     // debug = true;
0950
0951     EventOfCombCandidates &eoccs = m_event_of_comb_cands;
0952
0953     tbb::parallel_for_each(m_job->regions_begin(), m_job->regions_end(), [&](int region) {
0954       if (iteration_dir == SteeringParams::IT_BkwSearch && !m_job->steering_params(region).has_bksearch_plan()) {
0955         printf("No backward search plan for region %d\n", region);
0956         return;
0957       }
0958
0959       const RegionOfSeedIndices rosi(m_seedEtaSeparators, region);
0960
0961       // adaptive seeds per task based on the total estimated amount of work to divide among all threads
0962       const int adaptiveSPT = std::clamp(
0963           Config::numThreadsEvents * eoccs.size() / Config::numThreadsFinder + 1, 4, Config::numSeedsPerTask);
0964       dprint("adaptiveSPT " << adaptiveSPT << " fill " << rosi.count() << "/" << eoccs.size() << " region " << region);
0965
0966       tbb::parallel_for(rosi.tbb_blk_rng_std(adaptiveSPT), [&](const tbb::blocked_range<int> &seeds) {
0967         auto cloner = g_exe_ctx.m_cloners.makeOrGet();
0968         auto mkfndr = g_exe_ctx.m_finders.makeOrGet();
0969
0970         cloner->setup(m_job->params());
0971
0972         // loop over layers
0973         find_tracks_in_layers(*cloner, mkfndr.get(), iteration_dir, seeds.begin(), seeds.end(), region);
0974
0975         mkfndr->release();
0976         cloner->release();
0977       });
0978     });
0979
0980     // debug = false;
0981   }
0982
0983   void MkBuilder::find_tracks_in_layers(CandCloner &cloner,
0984                                         MkFinder *mkfndr,
0985                                         SteeringParams::IterationType_e iteration_dir,
0986                                         const int start_seed,
0987                                         const int end_seed,
0988                                         const int region) {
0989     EventOfCombCandidates &eoccs = m_event_of_comb_cands;
0990     const TrackerInfo &trk_info = m_job->m_trk_info;
0991     const SteeringParams &st_par = m_job->steering_params(region);
0992     const IterationParams &params = m_job->params();
0993     const PropagationConfig &prop_config = trk_info.prop_config();
0994
0995     const int n_seeds = end_seed - start_seed;
0996
0997     std::vector<std::pair<int, int>> seed_cand_idx;
0998     std::vector<UpdateIndices> seed_cand_update_idx, seed_cand_overlap_idx;
0999     seed_cand_idx.reserve(n_seeds * params.maxCandsPerSeed);
1000     seed_cand_update_idx.reserve(n_seeds * params.maxCandsPerSeed);
1001     seed_cand_overlap_idx.reserve(n_seeds * params.maxCandsPerSeed);
1002
1003     std::vector<std::vector<TrackCand>> extra_cands(n_seeds);
1004     for (int ii = 0; ii < n_seeds; ++ii)
1005       extra_cands[ii].reserve(params.maxCandsPerSeed);
1006
1007     cloner.begin_eta_bin(&eoccs, &seed_cand_update_idx, &seed_cand_overlap_idx, &extra_cands, start_seed, n_seeds);
1008
1009     // Loop over layers, starting from after the seed.
1010
1011     auto layer_plan_it = st_par.make_iterator(iteration_dir);
1012
1013     dprintf("Made iterator for %d, first layer=%d ... end layer=%d\n",
1014             iteration_dir,
1015             layer_plan_it.layer(),
1016             layer_plan_it.last_layer());
1017
1018     assert(layer_plan_it.is_pickup_only());
1019
1020     int curr_layer = layer_plan_it.layer(), prev_layer;
1021
1022     dprintf(
1023         "\nMkBuilder::find_tracks_in_layers region=%d, seed_pickup_layer=%d, first_layer=%d; start_seed=%d, "
1024         "end_seed=%d\n",
1025         region,
1026         curr_layer,
1027         layer_plan_it.next_layer(),
1028         start_seed,
1029         end_seed);
1030
1031     auto &iter_params = (iteration_dir == SteeringParams::IT_BkwSearch) ? m_job->m_iter_config.m_backward_params
1032                                                                         : m_job->m_iter_config.m_params;
1033
1034     // Loop over layers according to plan.
1035     while (++layer_plan_it) {
1036       prev_layer = curr_layer;
1037       curr_layer = layer_plan_it.layer();
1038       mkfndr->setup(prop_config,
1039                     m_job->m_iter_config,
1040                     iter_params,
1041                     m_job->m_iter_config.m_layer_configs[curr_layer],
1042                     st_par,
1043                     m_job->get_mask_for_layer(curr_layer),
1044                     m_event,
1045                     region,
1046                     m_job->m_in_fwd);
1047
1048       const bool pickup_only = layer_plan_it.is_pickup_only();
1049
1050       const LayerInfo &layer_info = trk_info.layer(curr_layer);
1051       const LayerOfHits &layer_of_hits = m_job->m_event_of_hits[curr_layer];
1052       const FindingFoos &fnd_foos = FindingFoos::get_finding_foos(layer_info.is_barrel());
1053
1054       dprintf("\n\n* Processing layer %d, %s\n\n", curr_layer, pickup_only ? "pickup only" : "full finding");
1055       mkfndr->begin_layer(layer_of_hits);
1056
1057       const int theEndCand = find_tracks_unroll_candidates(
1058           seed_cand_idx, start_seed, end_seed, curr_layer, prev_layer, pickup_only, iteration_dir);
1059
1060       dprintf("  Number of candidates to process: %d, nHits in layer: %d\n", theEndCand, layer_of_hits.nHits());
1061
1062       // Don't bother messing with the clone engine if there are no candidates
1063       // (actually it crashes, so this protection is needed).
1064       // If there are no cands on this iteration, there won't be any later on either,
1065       // by the construction of the seed_cand_idx vector.
1066       // XXXXMT There might be cases in endcap where all tracks will miss the
1067       // next layer, but only relevant if we do geometric selection before.
1068
1069       if (pickup_only || theEndCand == 0)
1070         continue;
1071
1072       cloner.begin_layer(curr_layer);
1073
1074       //vectorized loop
1075       for (int itrack = 0; itrack < theEndCand; itrack += NN) {
1076         const int end = std::min(itrack + NN, theEndCand);
1077
1078 #ifdef DEBUG
1079         dprintf("\nProcessing track=%d, start_seed=%d, n_seeds=%d, theEndCand=%d, end=%d, nn=%d, end_eq_tec=%d\n",
1080                 itrack,
1081                 start_seed,
1082                 n_seeds,
1083                 theEndCand,
1084                 end,
1085                 end - itrack,
1086                 end == theEndCand);
1087         dprintf("  (seed,cand): ");
1088         for (int i = itrack; i < end; ++i)
1089           dprintf("(%d,%d)  ", seed_cand_idx[i].first, seed_cand_idx[i].second);
1090         dprintf("\n");
1091 #endif
1092
1093         mkfndr->inputTracksAndHitIdx(eoccs.refCandidates(), seed_cand_idx, itrack, end, false);
1094
1095 #ifdef DEBUG
1096         for (int i = itrack; i < end; ++i)
1097           dprintf("  track %d, idx %d is from seed %d\n", i, i - itrack, mkfndr->m_Label(i - itrack, 0, 0));
1098 #endif
1099
1100         // propagate to current layer
1101         mkfndr->clearFailFlag();
1102         (mkfndr->*fnd_foos.m_propagate_foo)(
1103             layer_info.propagate_to(), end - itrack, prop_config.finding_inter_layer_pflags);
1104
1105         dprint("now get hit range");
1106
1107         if (alwaysUseHitSelectionV2 || iter_params.useHitSelectionV2)
1108           mkfndr->selectHitIndicesV2(layer_of_hits, end - itrack);
1109         else
1110           mkfndr->selectHitIndices(layer_of_hits, end - itrack);
1111
1112         find_tracks_handle_missed_layers(
1113             mkfndr, layer_info, extra_cands, seed_cand_idx, region, start_seed, itrack, end);
1114
1115         // copy_out the propagated track params, errors only.
1116         // Do not, keep cands at last valid hit until actual update,
1117         // this requires change to propagation flags used in MkFinder::updateWithLastHit()
1118         // from intra-layer to inter-layer.
1119         // mkfndr->copyOutParErr(eoccs.refCandidates_nc(), end - itrack, true);
1120
1121         // For prop-to-plane propagate from the last hit, not layer center.
1122         if constexpr (Config::usePropToPlane) {
1123           mkfndr->inputTracksAndHitIdx(eoccs.refCandidates(), seed_cand_idx, itrack, end, false);
1124         }
1125
1126         dprint("make new candidates");
1127         cloner.begin_iteration();
1128
1129         mkfndr->findCandidatesCloneEngine(layer_of_hits, cloner, start_seed, end - itrack, fnd_foos);
1130
1131         cloner.end_iteration();
1132       }  //end of vectorized loop
1133
1134       cloner.end_layer();
1135
1136       // Update loop of best candidates. CandCloner prepares the list of those
1137       // that need update (excluding all those with negative last hit index).
1138       // This is split into two sections - candidates without overlaps and with overlaps.
1139       // On CMS PU-50 the ratio of those is ~ 65 : 35 over all iterations.
1140       // Note, overlap recheck is only enabled for some iterations, e.g. pixelLess.
1141
1142       const int theEndUpdater = seed_cand_update_idx.size();
1143
1144       for (int itrack = 0; itrack < theEndUpdater; itrack += NN) {
1145         const int end = std::min(itrack + NN, theEndUpdater);
1146
1147         mkfndr->inputTracksAndHits(eoccs.refCandidates(), layer_of_hits, seed_cand_update_idx, itrack, end, true);
1148
1149         mkfndr->updateWithLoadedHit(end - itrack, layer_of_hits, fnd_foos);
1150
1151         // copy_out the updated track params, errors only (hit-idcs and chi2 already set)
1152         mkfndr->copyOutParErr(eoccs.refCandidates_nc(), end - itrack, false);
1153       }
1154
1155       const int theEndOverlapper = seed_cand_overlap_idx.size();
1156
1157       for (int itrack = 0; itrack < theEndOverlapper; itrack += NN) {
1158         const int end = std::min(itrack + NN, theEndOverlapper);
1159
1160         mkfndr->inputTracksAndHits(eoccs.refCandidates(), layer_of_hits, seed_cand_overlap_idx, itrack, end, true);
1161
1162         mkfndr->updateWithLoadedHit(end - itrack, layer_of_hits, fnd_foos);
1163
1164         mkfndr->copyOutParErr(eoccs.refCandidates_nc(), end - itrack, false);
1165
1166         mkfndr->inputOverlapHits(layer_of_hits, seed_cand_overlap_idx, itrack, end);
1167
1168         // XXXX Could also be calcChi2AndUpdate(), then copy-out would have to be done
1169         // below, choosing appropriate slot (with or without the overlap hit).
1170         // Probably in a dedicated MkFinder copyOutXyzz function.
1171         mkfndr->chi2OfLoadedHit(end - itrack, fnd_foos);
1172
1173         for (int ii = itrack; ii < end; ++ii) {
1174           const int fi = ii - itrack;
1175           TrackCand &tc = eoccs[seed_cand_overlap_idx[ii].seed_idx][seed_cand_overlap_idx[ii].cand_idx];
1176
1177           // XXXX For now we DO NOT use chi2 as this was how things were done before the post-update
1178           // chi2 check. To use it we should retune scoring function (might be even simpler).
1179           auto chi2Ovlp = mkfndr->m_Chi2[fi];
1180           if (mkfndr->m_FailFlag[fi] == 0 && chi2Ovlp >= 0.0f && chi2Ovlp <= 60.0f) {
1181             auto scoreCand =
1182                 getScoreCand(st_par.m_track_scorer, tc, true /*penalizeTailMissHits*/, true /*inFindCandidates*/);
1183             tc.addHitIdx(seed_cand_overlap_idx[ii].ovlp_idx, curr_layer, chi2Ovlp);
1184             tc.incOverlapCount();
1185             auto scoreCandOvlp = getScoreCand(st_par.m_track_scorer, tc, true, true);
1186             if (scoreCand > scoreCandOvlp)
1187               tc.popOverlap();
1188           }
1189         }
1190       }
1191
1192       // Check if cands are sorted, as expected.
1193 #ifdef DEBUG
1194       for (int iseed = start_seed; iseed < end_seed; ++iseed) {
1195         auto &cc = eoccs[iseed];
1196
1197         for (int i = 0; i < ((int)cc.size()) - 1; ++i) {
1198           if (cc[i].score() < cc[i + 1].score()) {
1199             printf("CloneEngine - NOT SORTED: layer=%d, iseed=%d (size=%lu)-- %d : %f smaller than %d : %f\n",
1200                    curr_layer,
1201                    iseed,
1202                    cc.size(),
1203                    i,
1204                    cc[i].score(),
1205                    i + 1,
1206                    cc[i + 1].score());
1207           }
1208         }
1209       }
1210 #endif
1211       mkfndr->end_layer();
1212     }  // end of layer loop
1213
1214     cloner.end_eta_bin();
1215
1216     // final sorting
1217     for (int iseed = start_seed; iseed < end_seed; ++iseed) {
1218       eoccs[iseed].mergeCandsAndBestShortOne(m_job->params(), st_par.m_track_scorer, true, true);
1219     }
1220   }
1221
1222   //==============================================================================
1223   // BackwardFit
1224   //==============================================================================
1225
1226 #ifdef DEBUG_FINAL_FIT
1227   namespace {
1228     // clang-format off
1229     void dprint_tcand(const TrackCand &t, int i) {
1230       dprintf("  %4d with q=%+d chi2=%7.3f pT=%7.3f eta=% 7.3f x=%.3f y=%.3f z=%.3f"
1231               " nHits=%2d  label=%4d findable=%d\n",
1232               i, t.charge(), t.chi2(), t.pT(), t.momEta(), t.x(), t.y(), t.z(),
1233               t.nFoundHits(), t.label(), t.isFindable());
1234       }
1235     // clang-format on
1236   }  // namespace
1237 #endif
1238
1239   void MkBuilder::backwardFitBH() {
1240     tbb::parallel_for_each(m_job->regions_begin(), m_job->regions_end(), [&](int region) {
1241       const RegionOfSeedIndices rosi(m_seedEtaSeparators, region);
1242
1243       tbb::parallel_for(rosi.tbb_blk_rng_vec(), [&](const tbb::blocked_range<int> &blk_rng) {
1244         auto mkfndr = g_exe_ctx.m_finders.makeOrGet();
1245
1246         RangeOfSeedIndices rng = rosi.seed_rng(blk_rng);
1247
1248         while (rng.valid()) {
1249           // final backward fit
1250           fit_cands_BH(mkfndr.get(), rng.m_beg, rng.m_end, region);
1251
1252           ++rng;
1253         }
1254       });
1255     });
1256   }
1257
1258   void MkBuilder::fit_cands_BH(MkFinder *mkfndr, int start_cand, int end_cand, int region) {
1259     const SteeringParams &st_par = m_job->steering_params(region);
1260     const PropagationConfig &prop_config = m_job->m_trk_info.prop_config();
1261     mkfndr->setup_bkfit(prop_config, st_par, m_event);
1262 #ifdef DEBUG_FINAL_FIT
1263     EventOfCombCandidates &eoccs = m_event_of_comb_cands;
1264     bool debug = true;
1265 #endif
1266
1267     for (int icand = start_cand; icand < end_cand; icand += NN) {
1268       const int end = std::min(icand + NN, end_cand);
1269
1270 #ifdef DEBUG_FINAL_FIT
1271       dprintf("Pre Final fit for %d - %d\n", icand, end);
1272       for (int i = icand; i < end; ++i) {
1273         dprint_tcand(eoccs[i][0], i);
1274       }
1275 #endif
1276
1277       bool chi_debug = false;
1278 #ifdef DEBUG_BACKWARD_FIT_BH
1279     redo_fit:
1280 #endif
1281
1282       // input candidate tracks
1283       mkfndr->bkFitInputTracks(m_tracks, icand, end);
1284
1285       // perform fit back to first layer on track
1286       mkfndr->bkFitFitTracksBH(m_job->m_event_of_hits, st_par, end - icand, chi_debug);
1287
1288       // now move one last time to PCA
1289       if (prop_config.backward_fit_to_pca) {
1290         mkfndr->bkFitPropTracksToPCA(end - icand);
1291       }
1292
1293 #ifdef DEBUG_BACKWARD_FIT_BH
1294       // Dump tracks with pT > 2 and chi2/dof > 20. Assumes MPT_SIZE=1.
1295       if (!chi_debug && 1.0f / mkfndr->m_Par[MkBase::iP].At(0, 3, 0) > 2.0f &&
1296           mkfndr->m_Chi2(0, 0, 0) / (eoccs[icand][0].nFoundHits() * 3 - 6) > 20.0f) {
1297         chi_debug = true;
1298 #ifdef MKFIT_STANDALONE
1299         printf("CHIHDR Event %d, Cand %3d, pT %f, chipdof %f ### NOTE x,y,z in cm, sigmas, deltas in mum ### !!!\n",
1300                m_event->evtID(),
1301 #else
1302         printf("CHIHDR Cand %3d, pT %f, chipdof %f ### NOTE x,y,z in cm, sigmas, deltas in mum ### !!!\n",
1303 #endif
1304                icand,
1305                1.0f / mkfndr->m_Par[MkBase::iP].At(0, 3, 0),
1306                mkfndr->m_Chi2(0, 0, 0) / (eoccs[icand][0].nFoundHits() * 3 - 6));
1307         // clang-format off
1308         printf("CHIHDR %3s %10s"
1309               " %10s %10s %10s %10s %11s %11s %11s"
1310               " %10s %10s %10s %10s %11s %11s %11s"
1311               " %10s %10s %10s %10s %10s %11s %11s\n",
1312               "lyr", "chi2",
1313               "x_h", "y_h", "z_h", "r_h", "sx_h", "sy_h", "sz_h",
1314               "x_t", "y_t", "z_t", "r_t", "sx_t", "sy_t", "sz_t",
1315               "pt", "phi", "theta", "phi_h", "phi_t", "d_xy", "d_z");
1316         // clang-format on
1317         goto redo_fit;
1318       }
1319 #endif
1320
1321       // copy out full set of info at last propagated position
1322       mkfndr->bkFitOutputTracks(m_tracks, icand, end, prop_config.backward_fit_to_pca);
1323
1324 #ifdef DEBUG_FINAL_FIT
1325       dprintf("Post Final fit for %d - %d\n", icand, end);
1326       for (int i = icand; i < end; ++i) {
1327         dprint_tcand(eoccs[i][0], i);
1328       }
1329 #endif
1330     }
1331   }
1332
1333   //------------------------------------------------------------------------------
1334
1335   void MkBuilder::backwardFit() {
1336     EventOfCombCandidates &eoccs = m_event_of_comb_cands;
1337
1338     tbb::parallel_for_each(m_job->regions_begin(), m_job->regions_end(), [&](int region) {
1339       const RegionOfSeedIndices rosi(m_seedEtaSeparators, region);
1340
1341       // adaptive seeds per task based on the total estimated amount of work to divide among all threads
1342       const int adaptiveSPT = std::clamp(
1343           Config::numThreadsEvents * eoccs.size() / Config::numThreadsFinder + 1, 4, Config::numSeedsPerTask);
1344       dprint("adaptiveSPT " << adaptiveSPT << " fill " << rosi.count() << "/" << eoccs.size() << " region " << region);
1345
1346       tbb::parallel_for(rosi.tbb_blk_rng_std(adaptiveSPT), [&](const tbb::blocked_range<int> &cands) {
1347         auto mkfndr = g_exe_ctx.m_finders.makeOrGet();
1348
1349         fit_cands(mkfndr.get(), cands.begin(), cands.end(), region);
1350       });
1351     });
1352   }
1353
1354   void MkBuilder::fit_cands(MkFinder *mkfndr, int start_cand, int end_cand, int region) {
1355     EventOfCombCandidates &eoccs = m_event_of_comb_cands;
1356     const SteeringParams &st_par = m_job->steering_params(region);
1357     const PropagationConfig &prop_config = m_job->m_trk_info.prop_config();
1358     mkfndr->setup_bkfit(prop_config, st_par, m_event);
1359
1360     int step = NN;
1361     for (int icand = start_cand; icand < end_cand; icand += step) {
1362       int end = std::min(icand + NN, end_cand);
1363
1364       bool chi_debug = false;
1365
1366 #ifdef DEBUG_FINAL_FIT
1367       bool debug = true;
1368       dprintf("Pre Final fit for %d - %d\n", icand, end);
1369       for (int i = icand; i < end; ++i) {
1370         dprint_tcand(eoccs[i][0], i);
1371       }
1372       chi_debug = true;
1373       static bool first = true;
1374       if (first) {
1375         // ./mkFit ... | perl -ne 'if (/^BKF_OVERLAP/) { s/^BKF_OVERLAP //og; print; }' > bkf_ovlp.rtt
1376         dprintf(
1377             "BKF_OVERLAP event/I:label/I:prod_type/I:is_findable/I:layer/I:is_stereo/I:is_barrel/I:"
1378             "pt/F:pt_cur/F:eta/F:phi/F:phi_cur/F:r_cur/F:z_cur/F:chi2/F:isnan/I:isfin/I:gtzero/I:hit_label/I:"
1379             "sx_t/F:sy_t/F:sz_t/F:d_xy/F:d_z/F\n");
1380         first = false;
1381       }
1382 #endif
1383
1384       // input tracks
1385       mkfndr->bkFitInputTracks(eoccs, icand, end);
1386
1387       // fit tracks back to first layer
1388       mkfndr->bkFitFitTracks(m_job->m_event_of_hits, st_par, end - icand, chi_debug);
1389
1390       // now move one last time to PCA
1391       if (prop_config.backward_fit_to_pca) {
1392         mkfndr->bkFitPropTracksToPCA(end - icand);
1393       }
1394
1395       mkfndr->bkFitOutputTracks(eoccs, icand, end, prop_config.backward_fit_to_pca);
1396
1397 #ifdef DEBUG_FINAL_FIT
1398       dprintf("Post Final fit for %d - %d\n", icand, end);
1399       for (int i = icand; i < end; ++i) {
1400         dprint_tcand(eoccs[i][0], i);
1401       }
1402 #endif
1403     }
1404     mkfndr->release();
1405   }
1406
1407 }  // end namespace mkfit