Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-27 23:38:50

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