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