MkFitCMS/standalone/buildtestMPlex.cc

0001 #include "RecoTracker/MkFitCMS/standalone/buildtestMPlex.h"
0002 #include "RecoTracker/MkFitCore/standalone/ConfigStandalone.h"
0003 #include "RecoTracker/MkFitCore/src/Matrix.h"
0004 #include "RecoTracker/MkFitCore/interface/MkBuilder.h"
0005 #include "RecoTracker/MkFitCMS/interface/MkStdSeqs.h"
0006 #include "RecoTracker/MkFitCMS/standalone/MkStandaloneSeqs.h"
0007
0008 #include "oneapi/tbb/parallel_for.h"
0009
0010 #include <memory>
0011
0012 namespace mkfit {
0013
0014   inline bool sortByHitsChi2(const std::pair<Track, TrackState> &cand1, const std::pair<Track, TrackState> &cand2) {
0015     if (cand1.first.nFoundHits() == cand2.first.nFoundHits())
0016       return cand1.first.chi2() < cand2.first.chi2();
0017
0018     return cand1.first.nFoundHits() > cand2.first.nFoundHits();
0019   }
0020
0021   inline bool sortByPhi(const Hit &hit1, const Hit &hit2) {
0022     return std::atan2(hit1.y(), hit1.x()) < std::atan2(hit2.y(), hit2.x());
0023   }
0024
0025   inline bool sortByEta(const Hit &hit1, const Hit &hit2) { return hit1.eta() < hit2.eta(); }
0026
0027   inline bool sortTracksByEta(const Track &track1, const Track &track2) { return track1.momEta() < track2.momEta(); }
0028
0029   inline bool sortTracksByPhi(const Track &track1, const Track &track2) { return track1.momPhi() < track2.momPhi(); }
0030
0031   struct sortTracksByPhiStruct {
0032     const std::vector<std::vector<Track>> &m_track_candidates;
0033
0034     sortTracksByPhiStruct(std::vector<std::vector<Track>> *track_candidates) : m_track_candidates(*track_candidates) {}
0035
0036     bool operator()(const std::pair<int, int> &track1, const std::pair<int, int> &track2) {
0037       return m_track_candidates[track1.first][track1.second].posPhi() <
0038              m_track_candidates[track2.first][track2.second].posPhi();
0039     }
0040   };
0041
0042   // within a layer with a "reasonable" geometry, ordering by Z is the same as eta
0043   inline bool sortByZ(const Hit &hit1, const Hit &hit2) { return hit1.z() < hit2.z(); }
0044
0045   //==============================================================================
0046   // NaN and Silly track parameter check
0047   //==============================================================================
0048
0049   namespace {
0050
0051     int check_nan_n_silly(TrackVec &tracks, const char *prefix) {
0052       int count = 0;
0053       for (auto &t : tracks) {
0054         if (t.hasSillyValues(Const::nan_n_silly_print_bad_cands_bkfit, false, prefix)) {
0055           ++count;
0056         }
0057       }
0058       return count;
0059     }
0060
0061     void check_nan_n_silly_candidates(Event &ev) {
0062       // MIMI -- nan_n_silly_per_layer_count is in MkBuilder, could be in MkJob.
0063       // if (Const::nan_n_silly_check_cands_every_layer)
0064       // {
0065       //   int sc = (int) ev.nan_n_silly_per_layer_count_;
0066       //   if (sc > 0)
0067       //     printf("Nan'n'Silly: Number of silly candidates over all layers = %d\n", sc);
0068       // }
0069       if (Const::nan_n_silly_check_cands_pre_bkfit) {
0070         int sc = check_nan_n_silly(ev.candidateTracks_, "Pre-bkfit silly check");
0071         if (sc > 0)
0072           printf("Nan'n'Silly: Number of silly pre-bkfit candidates = %d\n", sc);
0073       }
0074     }
0075
0076     void check_nan_n_silly_bkfit(Event &ev) {
0077       if (Const::nan_n_silly_check_cands_post_bkfit) {
0078         int sc = check_nan_n_silly(ev.fitTracks_, "Post-bkfit silly check");
0079         if (sc > 0)
0080           printf("Nan'n'Silly: Number of silly post-bkfit candidates = %d\n", sc);
0081       }
0082     }
0083
0084   }  // namespace
0085
0086   //==============================================================================
0087   // runBuildTestPlexDumbCMSSW
0088   //==============================================================================
0089
0090   void runBuildingTestPlexDumbCMSSW(Event &ev, const EventOfHits &eoh, MkBuilder &builder) {
0091     const IterationConfig &itconf = Config::ItrInfo[0];
0092
0093     MkJob job({Config::TrkInfo, itconf, eoh});
0094
0095     builder.begin_event(&job, &ev, __func__);
0096
0097     if (Config::sim_val_for_cmssw) {
0098       StdSeq::root_val_dumb_cmssw(&ev);
0099     }
0100
0101     builder.end_event();
0102   }
0103
0104   //==============================================================================
0105   // runBuildTestPlexBestHit
0106   //==============================================================================
0107
0108   double runBuildingTestPlexBestHit(Event &ev, const EventOfHits &eoh, MkBuilder &builder) {
0109     const IterationConfig &itconf = Config::ItrInfo[0];
0110
0111     const bool validation_on = (Config::sim_val || Config::quality_val);
0112
0113     if (validation_on) {
0114       TrackVec seeds1;
0115
0116       unsigned int algorithms[] = {4};  //only initialStep
0117
0118       for (auto const &s : ev.seedTracks_) {
0119         //keep seeds form the first iteration for processing
0120         if (std::find(algorithms, algorithms + 1, s.algoint()) != algorithms + 1)
0121           seeds1.push_back(s);
0122       }
0123       ev.seedTracks_.swap(seeds1);  //necessary for the validation - PrepareSeeds
0124       ev.relabel_bad_seedtracks();  //necessary for the validation - PrepareSeeds
0125     }
0126
0127     IterationMaskIfc mask_ifc;
0128
0129     // To disable hit-masks, pass nullptr in place of &mask_ifc to MkJob ctor
0130     // and optionally comment out ev.fill_hitmask_bool_vectors() call.
0131
0132     ev.fill_hitmask_bool_vectors(itconf.m_track_algorithm, mask_ifc.m_mask_vector);
0133
0134     MkJob job({Config::TrkInfo, itconf, eoh, &mask_ifc});
0135
0136     builder.begin_event(&job, &ev, __func__);
0137
0138     bool seeds_sorted = false;
0139     // CCCC builder.PrepareSeeds();
0140
0141     // EventOfCandidates event_of_cands;
0142     builder.find_tracks_load_seeds_BH(ev.seedTracks_, seeds_sorted);
0143
0144 #ifdef USE_VTUNE_PAUSE
0145     __SSC_MARK(0x111);  // use this to resume Intel SDE at the same point
0146     __itt_resume();
0147 #endif
0148
0149     double time = dtime();
0150
0151     builder.findTracksBestHit();
0152
0153     time = dtime() - time;
0154
0155 #ifdef USE_VTUNE_PAUSE
0156     __itt_pause();
0157     __SSC_MARK(0x222);  // use this to pause Intel SDE at the same point
0158 #endif
0159
0160     // Hack, get the tracks out.
0161     ev.candidateTracks_ = builder.ref_tracks();
0162
0163     // For best hit, the candidateTracks_ vector is the direct input to the backward fit so only need to do find_duplicates once
0164     if (Config::quality_val || Config::sim_val || Config::cmssw_val) {
0165       //Mark tracks as duplicates; if within CMSSW, remove duplicate tracks before backward fit
0166       if (Config::removeDuplicates) {
0167         StdSeq::find_duplicates(ev.candidateTracks_);
0168       }
0169     }
0170
0171     // now do backwards fit... do we want to time this section?
0172     if (Config::backwardFit) {
0173       builder.backwardFitBH();
0174       ev.fitTracks_ = builder.ref_tracks();
0175     }
0176
0177     if (Config::quality_val) {
0178       StdSeq::Quality qval;
0179       qval.quality_val(&ev);
0180     } else if (Config::sim_val || Config::cmssw_val) {
0181       StdSeq::root_val(&ev);
0182     }
0183
0184     builder.end_event();
0185
0186     // ev.print_tracks(ev.candidateTracks_, true);
0187
0188     return time;
0189   }
0190
0191   //==============================================================================
0192   // runBuildTestPlex Combinatorial: Standard TBB
0193   //==============================================================================
0194
0195   double runBuildingTestPlexStandard(Event &ev, const EventOfHits &eoh, MkBuilder &builder) {
0196     const IterationConfig &itconf = Config::ItrInfo[0];
0197
0198     const bool validation_on = (Config::sim_val || Config::quality_val);
0199
0200     if (validation_on) {
0201       TrackVec seeds1;
0202
0203       unsigned int algorithms[] = {4};  //only initialStep
0204
0205       for (auto const &s : ev.seedTracks_) {
0206         //keep seeds form the first iteration for processing
0207         if (std::find(algorithms, algorithms + 1, s.algoint()) != algorithms + 1)
0208           seeds1.push_back(s);
0209       }
0210       ev.seedTracks_.swap(seeds1);  //necessary for the validation - PrepareSeeds
0211       ev.relabel_bad_seedtracks();  //necessary for the validation - PrepareSeeds
0212     }
0213
0214     IterationMaskIfc mask_ifc;
0215
0216     // To disable hit-masks, pass nullptr in place of &mask_ifc to MkJob ctor
0217     // and optionally comment out ev.fill_hitmask_bool_vectors() call.
0218
0219     ev.fill_hitmask_bool_vectors(itconf.m_track_algorithm, mask_ifc.m_mask_vector);
0220
0221     MkJob job({Config::TrkInfo, itconf, eoh, &mask_ifc});
0222
0223     builder.begin_event(&job, &ev, __func__);
0224
0225     bool seeds_sorted = false;
0226     // CCCC builder.PrepareSeeds();
0227
0228     builder.find_tracks_load_seeds(ev.seedTracks_, seeds_sorted);
0229
0230 #ifdef USE_VTUNE_PAUSE
0231     __SSC_MARK(0x111);  // use this to resume Intel SDE at the same point
0232     __itt_resume();
0233 #endif
0234
0235     double time = dtime();
0236
0237     builder.findTracksStandard();
0238
0239     time = dtime() - time;
0240
0241 #ifdef USE_VTUNE_PAUSE
0242     __itt_pause();
0243     __SSC_MARK(0x222);  // use this to pause Intel SDE at the same point
0244 #endif
0245
0246     check_nan_n_silly_candidates(ev);
0247
0248     // first store candidate tracks
0249     builder.export_best_comb_cands(ev.candidateTracks_);
0250
0251     // now do backwards fit... do we want to time this section?
0252     if (Config::backwardFit) {
0253       // Using the TrackVec version until we home in on THE backward fit etc.
0254       // builder.backwardFit();
0255       builder.select_best_comb_cands();
0256       builder.backwardFitBH();
0257       ev.fitTracks_ = builder.ref_tracks();
0258
0259       check_nan_n_silly_bkfit(ev);
0260     }
0261
0262     StdSeq::handle_duplicates(&ev);
0263
0264     if (Config::quality_val) {
0265       StdSeq::Quality qval;
0266       qval.quality_val(&ev);
0267     } else if (Config::sim_val || Config::cmssw_val) {
0268       StdSeq::root_val(&ev);
0269     }
0270
0271     builder.end_event();
0272
0273     // ev.print_tracks(ev.candidateTracks_, true);
0274
0275     return time;
0276   }
0277
0278   //==============================================================================
0279   // runBuildTestPlex Combinatorial: CloneEngine TBB
0280   //==============================================================================
0281
0282   double runBuildingTestPlexCloneEngine(Event &ev, const EventOfHits &eoh, MkBuilder &builder) {
0283     const IterationConfig &itconf = Config::ItrInfo[0];
0284
0285     const bool validation_on = (Config::sim_val || Config::quality_val);
0286
0287     if (validation_on) {
0288       TrackVec seeds1;
0289
0290       unsigned int algorithms[] = {4};  //only initialStep
0291
0292       for (auto const &s : ev.seedTracks_) {
0293         //keep seeds form the first iteration for processing
0294         if (std::find(algorithms, algorithms + 1, s.algoint()) != algorithms + 1)
0295           seeds1.push_back(s);
0296       }
0297       ev.seedTracks_.swap(seeds1);  //necessary for the validation - PrepareSeeds
0298       ev.relabel_bad_seedtracks();  //necessary for the validation - PrepareSeeds
0299     }
0300
0301     IterationMaskIfc mask_ifc;
0302
0303     // To disable hit-masks, pass nullptr in place of &mask_ifc to MkJob ctor
0304     // and optionally comment out ev.fill_hitmask_bool_vectors() call.
0305
0306     ev.fill_hitmask_bool_vectors(itconf.m_track_algorithm, mask_ifc.m_mask_vector);
0307
0308     MkJob job({Config::TrkInfo, itconf, eoh, &mask_ifc});
0309
0310     builder.begin_event(&job, &ev, __func__);
0311
0312     bool seeds_sorted = false;
0313     // CCCC builder.PrepareSeeds();
0314
0315     builder.find_tracks_load_seeds(ev.seedTracks_, seeds_sorted);
0316
0317 #ifdef USE_VTUNE_PAUSE
0318     __SSC_MARK(0x111);  // use this to resume Intel SDE at the same point
0319     __itt_resume();
0320 #endif
0321
0322     double time = dtime();
0323
0324     builder.findTracksCloneEngine();
0325
0326     time = dtime() - time;
0327
0328 #ifdef USE_VTUNE_PAUSE
0329     __itt_pause();
0330     __SSC_MARK(0x222);  // use this to pause Intel SDE at the same point
0331 #endif
0332
0333     check_nan_n_silly_candidates(ev);
0334
0335     // first store candidate tracks - needed for BH backward fit and root_validation
0336     builder.export_best_comb_cands(ev.candidateTracks_);
0337
0338     // now do backwards fit... do we want to time this section?
0339     if (Config::backwardFit) {
0340       // a) TrackVec version:
0341       builder.select_best_comb_cands();
0342       builder.backwardFitBH();
0343       ev.fitTracks_ = builder.ref_tracks();
0344
0345       // b) Version that runs on CombCand / TrackCand
0346       // builder.backwardFit();
0347       // builder.quality_store_tracks(ev.fitTracks_);
0348
0349       check_nan_n_silly_bkfit(ev);
0350     }
0351
0352     StdSeq::handle_duplicates(&ev);
0353
0354     // validation section
0355     if (Config::quality_val) {
0356       StdSeq::Quality qval;
0357       qval.quality_val(&ev);
0358     } else if (Config::sim_val || Config::cmssw_val) {
0359       StdSeq::root_val(&ev);
0360     }
0361
0362     builder.end_event();
0363
0364     // ev.print_tracks(ev.candidateTracks_, true);
0365
0366     return time;
0367   }
0368
0369   //==============================================================================
0370   // runBtpCe_MultiIter
0371   //
0372   // Prototype for running multiple iterations, sequentially, using the same builder.
0373   // For cmmsw seeds
0374   //
0375   // There is, in general, a mess in how tracks are processed, marked, or copied out
0376   // in various validation scenarios and export flags.
0377   //
0378   // In particular, MkBuilder::PrepareSeeds does a lot of things to whole / complete
0379   // event,seedTracks_ -- probably this would need to be split into common / and
0380   // per-iteration part.
0381   // - MkBuilder::prep_*** functions also mostly do not belong there (prep_sim is called from
0382   //   PrepareSeeds() for cmssw seeds).
0383   //
0384   // At this point we need to think about what should happen to Event before all the iterations and
0385   // after all the iterations ... from the Validation perspective.
0386   // And if we care about doing too muich work for seeds that will never get processed.
0387   //==============================================================================
0388
0389   std::vector<double> runBtpCe_MultiIter(Event &ev, const EventOfHits &eoh, MkBuilder &builder, int n) {
0390     std::vector<double> timevec;
0391     if (n <= 0)
0392       return timevec;
0393     timevec.resize(n + 1, 0.0);
0394
0395     const bool validation_on = (Config::sim_val || Config::quality_val);
0396
0397     TrackVec seeds_used;
0398     TrackVec seeds1;
0399
0400     unsigned int algorithms[] = {4, 22, 23, 5, 24, 7, 8, 9, 10, 6};  //9 iterations
0401
0402     if (validation_on) {
0403       for (auto const &s : ev.seedTracks_) {
0404         //keep seeds form the first n iterations for processing
0405         if (std::find(algorithms, algorithms + n, s.algoint()) != algorithms + n)
0406           seeds1.push_back(s);
0407       }
0408       ev.seedTracks_.swap(seeds1);  //necessary for the validation - PrepareSeeds
0409       ev.relabel_bad_seedtracks();  //necessary for the validation - PrepareSeeds
0410     }
0411
0412     IterationMaskIfc mask_ifc;
0413     TrackVec seeds;
0414     TrackVec tmp_tvec;
0415
0416     for (int it = 0; it <= n - 1; ++it) {
0417       const IterationConfig &itconf = Config::ItrInfo[it];
0418
0419       // To disable hit-masks, pass nullptr in place of &mask_ifc to MkJob ctor
0420       // and optionally comment out ev.fill_hitmask_bool_vectors() call.
0421
0422       ev.fill_hitmask_bool_vectors(itconf.m_track_algorithm, mask_ifc.m_mask_vector);
0423
0424       MkJob job({Config::TrkInfo, itconf, eoh, &mask_ifc});
0425
0426       builder.begin_event(&job, &ev, __func__);
0427
0428       {  // We could partition seeds once, store beg, end for each iteration in a map or vector.
0429         seeds.clear();
0430         int nc = 0;
0431         for (auto &s : ev.seedTracks_) {
0432           if (s.algoint() == itconf.m_track_algorithm) {
0433             if (itconf.m_requires_seed_hit_sorting) {
0434               s.sortHitsByLayer();
0435             }
0436             seeds.push_back(s);
0437             ++nc;
0438           } else if (nc > 0)
0439             break;
0440         }
0441       }
0442
0443       if (itconf.m_requires_dupclean_tight)
0444         StdSeq::clean_cms_seedtracks_iter(&seeds, itconf, eoh.refBeamSpot());
0445
0446       builder.seed_post_cleaning(seeds);
0447
0448       // Add protection in case no seeds are found for iteration
0449       if (seeds.size() <= 0)
0450         continue;
0451
0452       builder.find_tracks_load_seeds(seeds, itconf.m_requires_dupclean_tight);
0453
0454       double time = dtime();
0455
0456       builder.findTracksCloneEngine();
0457
0458       timevec[it] = dtime() - time;
0459       timevec[n] += timevec[it];
0460
0461       // Print min and max size of hots vectors of CombCands.
0462       // builder.find_min_max_hots_size();
0463
0464       if (validation_on)
0465         seeds_used.insert(seeds_used.end(), seeds.begin(), seeds.end());  //cleaned seeds need to be stored somehow
0466
0467       using Algo = TrackBase::TrackAlgorithm;
0468       if (itconf.m_requires_quality_filter && Algo(itconf.m_track_algorithm) != Algo::detachedTripletStep) {
0469         if (Algo(itconf.m_track_algorithm) == Algo::pixelPairStep) {
0470           builder.filter_comb_cands([&](const TrackCand &t) { return StdSeq::qfilter_n_hits_pixseed(t, 3); });
0471         } else if (Algo(itconf.m_track_algorithm) == Algo::pixelLessStep) {
0472           builder.filter_comb_cands(
0473               [&](const TrackCand &t) { return StdSeq::qfilter_pixelLessFwd(t, eoh.refBeamSpot(), Config::TrkInfo); });
0474         } else {
0475           builder.filter_comb_cands(
0476               [&](const TrackCand &t) { return StdSeq::qfilter_n_hits(t, itconf.m_params.minHitsQF); });
0477         }
0478       }
0479
0480       builder.select_best_comb_cands();
0481
0482       {
0483         builder.export_tracks(tmp_tvec);
0484         StdSeq::find_and_remove_duplicates(tmp_tvec, itconf);
0485         ev.candidateTracks_.reserve(ev.candidateTracks_.size() + tmp_tvec.size());
0486         for (auto &&t : tmp_tvec)
0487           ev.candidateTracks_.emplace_back(std::move(t));
0488         tmp_tvec.clear();
0489       }
0490
0491       // now do backwards fit... do we want to time this section?
0492       if (Config::backwardFit) {
0493         // a) TrackVec version:
0494         // builder.backwardFitBH();
0495
0496         // b) Version that runs on CombCand / TrackCand
0497         const bool do_backward_search = Config::backwardSearch && itconf.m_backward_search;
0498
0499         // We copy seed-hits into Candidates ... now we have to remove them so backward fit stops
0500         // before reaching seeding region. Ideally, we wouldn't add them in the first place but
0501         // if we want to export full tracks above we need to hold on to them (alternatively, we could
0502         // have a pointer to seed track in CombCandidate and copy them from there).
0503         if (do_backward_search) {
0504           builder.compactifyHitStorageForBestCand(itconf.m_backward_drop_seed_hits, itconf.m_backward_fit_min_hits);
0505         }
0506
0507         builder.backwardFit();
0508
0509         if (do_backward_search) {
0510           builder.beginBkwSearch();
0511           builder.findTracksCloneEngine(SteeringParams::IT_BkwSearch);
0512           builder.endBkwSearch();
0513         }
0514
0515         if (itconf.m_requires_quality_filter && (Algo(itconf.m_track_algorithm) == Algo::detachedTripletStep ||
0516                                                  Algo(itconf.m_track_algorithm) == Algo::pixelLessStep)) {
0517           if (Algo(itconf.m_track_algorithm) == Algo::detachedTripletStep) {
0518             builder.filter_comb_cands(
0519                 [&](const TrackCand &t) { return StdSeq::qfilter_n_layers(t, eoh.refBeamSpot(), Config::TrkInfo); });
0520           } else if (Algo(itconf.m_track_algorithm) == Algo::pixelLessStep) {
0521             builder.filter_comb_cands([&](const TrackCand &t) {
0522               return StdSeq::qfilter_pixelLessBkwd(t, eoh.refBeamSpot(), Config::TrkInfo);
0523             });
0524           }
0525         }
0526
0527         builder.filter_comb_cands([&](const TrackCand &t) { return StdSeq::qfilter_nan_n_silly(t); });
0528
0529         builder.select_best_comb_cands(true);  // true -> clear m_tracks as they were already filled once above
0530
0531         StdSeq::find_and_remove_duplicates(builder.ref_tracks_nc(), itconf);
0532         builder.export_tracks(ev.fitTracks_);
0533       }
0534
0535       builder.end_event();
0536     }
0537
0538     // MIMI - Fake back event pointer for final processing (that should be done elsewhere)
0539     MkJob job({Config::TrkInfo, Config::ItrInfo[0], eoh});
0540     builder.begin_event(&job, &ev, __func__);
0541
0542     if (validation_on) {
0543       StdSeq::prep_simtracks(&ev);
0544       //swap for the cleaned seeds
0545       ev.seedTracks_.swap(seeds_used);
0546     }
0547
0548     check_nan_n_silly_candidates(ev);
0549
0550     if (Config::backwardFit)
0551       check_nan_n_silly_bkfit(ev);
0552
0553     // validation section
0554     if (Config::quality_val) {
0555       StdSeq::Quality qval;
0556       qval.quality_val(&ev);
0557     } else if (Config::sim_val || Config::cmssw_val) {
0558       StdSeq::root_val(&ev);
0559     }
0560
0561     // ev.print_tracks(ev.candidateTracks_, true);
0562
0563     // MIMI Unfake.
0564     builder.end_event();
0565
0566     // In CMSSW runOneIter we now release memory for comb-cands:
0567     builder.release_memory();
0568
0569     return timevec;
0570   }
0571
0572 }  // end namespace mkfit