Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-01-23 02:42:46

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