Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 11:22:29

0001 #include "RecoTracker/MkFitCMS/standalone/Shell.h"
0002 
0003 #include "RecoTracker/MkFitCore/src/Debug.h"
0004 
0005 // #include "RecoTracker/MkFitCore/src/Matriplex/MatriplexCommon.h"
0006 
0007 #include "RecoTracker/MkFitCMS/interface/runFunctions.h"
0008 
0009 #include "RecoTracker/MkFitCore/interface/HitStructures.h"
0010 #include "RecoTracker/MkFitCore/interface/MkBuilder.h"
0011 #include "RecoTracker/MkFitCore/src/MkFitter.h"
0012 #include "RecoTracker/MkFitCMS/interface/MkStdSeqs.h"
0013 #include "RecoTracker/MkFitCMS/standalone/MkStandaloneSeqs.h"
0014 
0015 #include "RecoTracker/MkFitCore/interface/Config.h"
0016 #include "RecoTracker/MkFitCore/standalone/ConfigStandalone.h"
0017 
0018 #include "RecoTracker/MkFitCore/standalone/Event.h"
0019 
0020 #ifndef NO_ROOT
0021 #include "TROOT.h"
0022 #include "TRint.h"
0023 #endif
0024 
0025 #include "oneapi/tbb/task_arena.h"
0026 
0027 #include <vector>
0028 
0029 // clang-format off
0030 
0031 namespace {
0032   constexpr int algos[] = {4, 22, 23, 5, 24, 7, 8, 9, 10, 6};  // 10 iterations
0033   constexpr int n_algos = sizeof(algos) / sizeof(int);
0034 
0035   const char* b2a(bool b) { return b ? "true" : "false"; }
0036 }
0037 
0038 namespace mkfit {
0039 
0040   Shell::Shell(std::vector<DeadVec> &dv, const std::string &in_file, int start_ev)
0041     : m_deadvectors(dv)
0042   {
0043     m_eoh = new EventOfHits(Config::TrkInfo);
0044     m_builder = new MkBuilder(Config::silent);
0045 
0046     m_backward_fit = Config::backwardFit;
0047 
0048     m_data_file = new DataFile;
0049     m_evs_in_file = m_data_file->openRead(in_file, Config::TrkInfo.n_layers());
0050 
0051     m_event = new Event(0, Config::TrkInfo.n_layers());
0052     GoToEvent(start_ev);
0053   }
0054 
0055   void Shell::Run() {
0056 #ifndef NO_ROOT
0057     std::vector<const char *> argv = { "mkFit", "-l" };
0058     int argc = argv.size();
0059     TRint rint("mkFit-shell", &argc, const_cast<char**>(argv.data()));
0060 
0061     char buf[256];
0062     sprintf(buf, "mkfit::Shell &s = * (mkfit::Shell*) %p;", this);
0063     gROOT->ProcessLine(buf);
0064     printf("Shell &s variable is set\n");
0065 
0066     rint.Run(true);
0067     printf("Shell::Run finished\n");
0068 #else
0069     printf("Shell::Run() no root, we rot -- erroring out. Recompile with WITH_ROOT=1\n");
0070     exit(1);
0071 #endif
0072   }
0073 
0074   void Shell::Status() {
0075     printf("On event %d, selected iteration index %d, algo %d - %s\n"
0076           "  debug = %s, use_dead_modules = %s\n"
0077            "  clean_seeds = %s, backward_fit = %s, remove_duplicates = %s\n",
0078            m_event->evtID(), m_it_index, algos[m_it_index], TrackBase::algoint_to_cstr(algos[m_it_index]),
0079            b2a(g_debug), b2a(Config::useDeadModules),
0080            b2a(m_clean_seeds), b2a(m_backward_fit), b2a(m_remove_duplicates));
0081   }
0082 
0083   //===========================================================================
0084   // Event navigation / processing
0085   //===========================================================================
0086 
0087   void Shell::GoToEvent(int eid) {
0088     if (eid < 1) {
0089       fprintf(stderr, "Requested event %d is less than 1 -- 1 is the first event, %d is total number of events in file\n",
0090              eid, m_evs_in_file);
0091       throw std::runtime_error("event out of range");
0092     }
0093     if (eid > m_evs_in_file) {
0094       fprintf(stderr, "Requested event %d is grater than total number of events in file %d\n",
0095              eid, m_evs_in_file);
0096       throw std::runtime_error("event out of range");
0097     }
0098 
0099     int pos = m_event->evtID();
0100     if (eid > pos) {
0101       m_data_file->skipNEvents(eid - pos - 1);
0102     } else {
0103       m_data_file->rewind();
0104       m_data_file->skipNEvents(eid - 1);
0105     }
0106     m_event->reset(eid);
0107     m_event->read_in(*m_data_file);
0108     StdSeq::loadHitsAndBeamSpot(*m_event, *m_eoh);
0109     if (Config::useDeadModules) {
0110       StdSeq::loadDeads(*m_eoh, m_deadvectors);
0111     }
0112 
0113     printf("At event %d\n", eid);
0114   }
0115 
0116   void Shell::NextEvent(int skip) {
0117     GoToEvent(m_event->evtID() + skip);
0118   }
0119 
0120   void Shell::ProcessEvent(SeedSelect_e seed_select, int selected_seed, int count) {
0121     // count is only used for SS_IndexPreCleaning and SS_IndexPostCleaning.
0122     //       There are no checks for upper bounds, ie, if requested seeds beyond the first one exist.
0123 
0124     const IterationConfig &itconf = Config::ItrInfo[m_it_index];
0125     IterationMaskIfc mask_ifc;
0126     m_event->fill_hitmask_bool_vectors(itconf.m_track_algorithm, mask_ifc.m_mask_vector);
0127 
0128     m_seeds.clear();
0129     m_tracks.clear();
0130 
0131     {
0132       int n_algo = 0; // seeds are grouped by algo
0133       for (auto &s : m_event->seedTracks_) {
0134         if (s.algoint() == itconf.m_track_algorithm) {
0135           if (seed_select == SS_UseAll || seed_select == SS_IndexPostCleaning) {
0136             m_seeds.push_back(s);
0137           } else if (seed_select == SS_Label && s.label() == selected_seed) {
0138             m_seeds.push_back(s);
0139             break;
0140           } else if (seed_select == SS_IndexPreCleaning && n_algo >= selected_seed) {
0141             m_seeds.push_back(s);
0142             if (--count <= 0)
0143               break;
0144           }
0145           ++n_algo;
0146         } else if (n_algo > 0)
0147           break;
0148       }
0149     }
0150 
0151     printf("Shell::ProcessEvent running over %d seeds\n", (int) m_seeds.size());
0152 
0153     // Equivalent to run_OneIteration(...) without MkBuilder::release_memory().
0154     // If seed_select == SS_IndexPostCleaning the given seed is picked after cleaning.
0155     {
0156       const TrackerInfo &trackerInfo = Config::TrkInfo;
0157       const EventOfHits &eoh = *m_eoh;
0158       const IterationMaskIfcBase &it_mask_ifc = mask_ifc;
0159       MkBuilder &builder = *m_builder;
0160       TrackVec &seeds = m_seeds;
0161       TrackVec &out_tracks = m_tracks;
0162       bool do_seed_clean = m_clean_seeds;
0163       bool do_backward_fit = m_backward_fit;
0164       bool do_remove_duplicates = m_remove_duplicates;
0165 
0166       MkJob job({trackerInfo, itconf, eoh, eoh.refBeamSpot(), &it_mask_ifc});
0167 
0168       builder.begin_event(&job, m_event, __func__);
0169 
0170       // Seed cleaning not done on all iterations.
0171       do_seed_clean = m_clean_seeds && itconf.m_seed_cleaner;
0172 
0173       if (do_seed_clean)
0174         itconf.m_seed_cleaner(seeds, itconf, eoh.refBeamSpot());
0175 
0176       // Check nans in seeds -- this should not be needed when Slava fixes
0177       // the track parameter coordinate transformation.
0178       builder.seed_post_cleaning(seeds);
0179 
0180       if (seed_select == SS_IndexPostCleaning) {
0181         if (selected_seed >= 0 && selected_seed < (int)seeds.size()) {
0182           for (int i = 0; i < count; ++i)
0183             seeds[i] = seeds[selected_seed + i];
0184           seeds.resize(count);
0185         } else {
0186           seeds.clear();
0187         }
0188       }
0189 
0190       if (seeds.empty()) {
0191         if (seed_select != SS_UseAll)
0192           printf("Shell::ProcessEvent requested seed not found among seeds of the selected iteration.\n");
0193         else
0194           printf("Shell::ProcessEvent no seeds found.\n");
0195         return;
0196       }
0197 
0198       if (itconf.m_requires_seed_hit_sorting) {
0199         for (auto &s : seeds)
0200           s.sortHitsByLayer();  // sort seed hits for the matched hits (I hope it works here)
0201       }
0202 
0203       builder.find_tracks_load_seeds(seeds, do_seed_clean);
0204 
0205       builder.findTracksCloneEngine();
0206 
0207       printf("Shell::ProcessEvent post fwd search: %d comb-cands\n", builder.ref_eocc().size());
0208 
0209       // Pre backward-fit filtering.
0210       filter_candidates_func pre_filter;
0211       if (do_backward_fit && itconf.m_pre_bkfit_filter)
0212         pre_filter = [&](const TrackCand &tc, const MkJob &jb) -> bool {
0213           return itconf.m_pre_bkfit_filter(tc, jb) && StdSeq::qfilter_nan_n_silly<TrackCand>(tc, jb);
0214         };
0215       else if (itconf.m_pre_bkfit_filter)
0216         pre_filter = itconf.m_pre_bkfit_filter;
0217       else if (do_backward_fit)
0218         pre_filter = StdSeq::qfilter_nan_n_silly<TrackCand>;
0219       // pre_filter can be null if we are not doing backward fit as nan_n_silly will be run below.
0220       if (pre_filter)
0221         builder.filter_comb_cands(pre_filter, true);
0222 
0223       printf("Shell::ProcessEvent post pre-bkf-filter (%s) and nan-filter (%s) filter: %d comb-cands\n",
0224              b2a(bool(itconf.m_pre_bkfit_filter)), b2a(do_backward_fit), builder.ref_eocc().size());
0225 
0226       job.switch_to_backward();
0227 
0228       if (do_backward_fit) {
0229         if (itconf.m_backward_search) {
0230           builder.compactifyHitStorageForBestCand(itconf.m_backward_drop_seed_hits, itconf.m_backward_fit_min_hits);
0231         }
0232 
0233         builder.backwardFit();
0234 
0235         if (itconf.m_backward_search) {
0236           builder.beginBkwSearch();
0237           builder.findTracksCloneEngine(SteeringParams::IT_BkwSearch);
0238         }
0239 
0240         printf("Shell::ProcessEvent post backward fit / search: %d comb-cands\n", builder.ref_eocc().size());
0241       }
0242 
0243       // Post backward-fit filtering.
0244       filter_candidates_func post_filter;
0245       if (do_backward_fit && itconf.m_post_bkfit_filter)
0246         post_filter = [&](const TrackCand &tc, const MkJob &jb) -> bool {
0247           return itconf.m_post_bkfit_filter(tc, jb) && StdSeq::qfilter_nan_n_silly<TrackCand>(tc, jb);
0248         };
0249       else
0250         post_filter = StdSeq::qfilter_nan_n_silly<TrackCand>;
0251       // post_filter is always at least doing nan_n_silly filter.
0252       builder.filter_comb_cands(post_filter, true);
0253 
0254       printf("Shell::ProcessEvent post post-bkf-filter (%s) and nan-filter (true): %d comb-cands\n",
0255              b2a(do_backward_fit && itconf.m_post_bkfit_filter), builder.ref_eocc().size());
0256 
0257       if (do_backward_fit && itconf.m_backward_search)
0258         builder.endBkwSearch();
0259 
0260       builder.export_best_comb_cands(out_tracks, true);
0261 
0262       if (do_remove_duplicates && itconf.m_duplicate_cleaner) {
0263         itconf.m_duplicate_cleaner(out_tracks, itconf);
0264       }
0265 
0266       printf("Shell::ProcessEvent post remove-duplicates: %d comb-cands\n", (int) out_tracks.size());
0267 
0268       builder.end_event();
0269     }
0270 
0271     printf("Shell::ProcessEvent found %d tracks, number of seeds at end %d\n",
0272            (int) m_tracks.size(), (int) m_seeds.size());
0273   }
0274 
0275   //===========================================================================
0276   // Iteration selection
0277   //===========================================================================
0278 
0279   void Shell::SelectIterationIndex(int itidx) {
0280     if (itidx < 0 || itidx >= n_algos) {
0281       fprintf(stderr, "Requested iteration index out of range [%d, %d)", 0, n_algos);
0282       throw std::runtime_error("iteration index out of range");
0283     }
0284     m_it_index = itidx;
0285   }
0286 
0287   void Shell::SelectIterationAlgo(int algo) {
0288     for (int i = 0; i < n_algos; ++i) {
0289       if (algos[i] == algo) {
0290         m_it_index = i;
0291         return;
0292       }
0293     }
0294     fprintf(stderr, "Requested algo %d not found", algo);
0295     throw std::runtime_error("algo not found");
0296   }
0297 
0298   void Shell::PrintIterations() {
0299     printf("Shell::PrintIterations selected index = %d, %d iterations hardcoded as\n",
0300             m_it_index, n_algos);
0301     for (int i = 0; i < n_algos; ++i)
0302       printf("%d %2d %s\n", i, algos[i], TrackBase::algoint_to_cstr(algos[i]));
0303   }
0304 
0305   //===========================================================================
0306   // Flags / status setters
0307   //===========================================================================
0308 
0309   void Shell::SetDebug(bool b) { g_debug = b; }
0310   void Shell::SetCleanSeeds(bool b) { m_clean_seeds = b; }
0311   void Shell::SetBackwardFit(bool b) { m_backward_fit = b; }
0312   void Shell::SetRemoveDuplicates(bool b) { m_remove_duplicates = b; }
0313   void Shell::SetUseDeadModules(bool b) { Config::useDeadModules = b; }
0314 
0315   //===========================================================================
0316   // Analysis helpers
0317   //===========================================================================
0318 
0319   /*
0320     sim tracks are written to .bin files with a label equal to its own index.
0321     reco tracks labels are seed indices.
0322     seed labels are sim track indices
0323     --
0324     mkfit labels are seed indices in given iteration after cleaning (at seed load-time)
0325   */
0326 
0327   int Shell::LabelFromHits(Track &t, bool replace, float good_frac) {
0328     std::map<int, int> lab_cnt;
0329     for (int hi = 0; hi < t.nTotalHits(); ++hi) {
0330       auto hot = t.getHitOnTrack(hi);
0331       if (hot.index < 0)
0332         continue;
0333       const Hit &h = m_event->layerHits_[hot.layer][hot.index];
0334       int hl = m_event->simHitsInfo_[h.mcHitID()].mcTrackID_;
0335       if (hl >= 0)
0336         ++lab_cnt[hl];
0337     }
0338     int max_c = -1, max_l = -1;
0339     for (auto& x : lab_cnt) {
0340       if (x.second > max_c) {
0341         max_l = x.first;
0342         max_c = x.second;
0343       }
0344     }
0345     bool success = max_c >= good_frac * t.nFoundHits();
0346     int relabel = success ? max_l : -1;
0347     // printf("found_hits=%d, best_lab %d (%d hits), existing label=%d (replace flag=%s)\n",
0348     //        t.nFoundHits(), max_l, max_c, t.label(), b2a(replace));
0349     if (replace)
0350         t.setLabel(relabel);
0351     return relabel;
0352   }
0353 
0354   void Shell::FillByLabelMaps_CkfBase() {
0355     Event &ev = *m_event;
0356     const int track_algo = Config::ItrInfo[m_it_index].m_track_algorithm;
0357 
0358     m_ckf_map.clear();
0359     m_sim_map.clear();
0360     m_seed_map.clear();
0361     m_mkf_map.clear();
0362 
0363     // Pick ckf tracks with right algo and a good label.
0364     int rec_algo_match = 0;
0365     for (auto &t : ev.cmsswTracks_) {
0366       if (t.algoint() != track_algo)
0367         continue;
0368       ++rec_algo_match;
0369       int label = LabelFromHits(t, false, 0.5);
0370       if (label >= 0) {
0371         m_ckf_map.insert(std::make_pair(label, &t));
0372       }
0373     }
0374 
0375     // Pick sim tracks with labels found by ckf.
0376     for (auto &t : ev.simTracks_) {
0377       if (t.label() >= 0 && m_ckf_map.find(t.label()) != m_ckf_map.end()) {
0378         m_sim_map.insert(std::make_pair(t.label(), &t));
0379       }
0380     }
0381 
0382     // Pick seeds with right algo and a label found by ckf.
0383     for (auto &t : ev.seedTracks_) {
0384       if (t.algoint() == track_algo && t.label() >= 0 && m_ckf_map.find(t.label()) != m_ckf_map.end()) {
0385         m_seed_map.insert(std::make_pair(t.label(), &t));
0386       }
0387     }
0388     // Some seeds seem to be labeled -1, try fixing when not otherwise found.
0389     for (auto &t : ev.seedTracks_) {
0390       if (t.algoint() == track_algo && t.label() == -1) {
0391         int lab = LabelFromHits(t, true, 0.5);
0392         if (lab >= 0 && m_seed_map.find(lab) == m_seed_map.end()) {
0393           if (m_ckf_map.find(lab) != m_ckf_map.end()) {
0394             m_seed_map.insert(std::make_pair(t.label(), &t));
0395             printf("Saved seed with label -1 -> %d\n", lab);
0396           }
0397         }
0398       }
0399     }
0400 
0401     // Pick mkfit tracks, label by 
0402     for (auto &t : m_tracks) {
0403       int label = LabelFromHits(t, false, 0.5);
0404       if (label >= 0) {
0405         m_mkf_map.insert(std::make_pair(label, &t));
0406       }
0407     }
0408 
0409     printf("Shell::FillByLabelMaps reporting tracks with label >= 0, algo=%d (%s): "
0410            "ckf: %d of %d (same algo=%d)), sim: %d of %d, seed: %d of %d, mkfit: %d w/label of %d\n",
0411            track_algo, TrackBase::algoint_to_cstr(track_algo),
0412            (int) m_ckf_map.size(), (int) ev.cmsswTracks_.size(), rec_algo_match,
0413            (int) m_sim_map.size(), (int) ev.simTracks_.size(),
0414            (int) m_seed_map.size(), (int) m_seeds.size(),
0415            (int) m_mkf_map.size(), (int) m_tracks.size()
0416     );
0417   }
0418 
0419   bool Shell::CheckMkFitLayerPlanVsReferenceHits(const Track &mkft, const Track &reft, const std::string &name) {
0420     // Check if all hit-layers of a reference track reft are in the mkfit layer plan.
0421     // Returns true if all layers are in the plan.
0422     // String name is printed in front of label, expected to be SIMK or CKF.
0423     const IterationConfig &itconf = Config::ItrInfo[m_it_index];
0424     auto lp = itconf.m_steering_params[ mkft.getEtaRegion() ].m_layer_plan;
0425     bool ret = true;
0426     for (int hi = 0; hi < reft.nTotalHits(); ++hi) {
0427       auto hot = reft.getHitOnTrack(hi);
0428       if (std::find_if(lp.begin(), lp.end(), [=](auto &x){ return x.m_layer == hot.layer; }) == lp.end())
0429       {
0430         printf("CheckMkfLayerPlanVsCkfHits: layer %d not in layer plan for region %d, %s label=%d\n",
0431                 hot.layer, mkft.getEtaRegion(), name.c_str(), reft.label());
0432         ret = false;
0433       }
0434     }
0435     return ret;
0436   }
0437 
0438   //===========================================================================
0439   // Analysis drivers / main functions / Comparators
0440   //===========================================================================
0441 
0442   void Shell::Compare() {
0443     Event &ev = *m_event;
0444     const IterationConfig &itconf = Config::ItrInfo[m_it_index];
0445 
0446     FillByLabelMaps_CkfBase();
0447 
0448     printf("------------------------------------------------------\n");
0449 
0450     const bool print_all_def = false;
0451     int mkf_cnt=0, less_hits=0, more_hits=0;
0452 
0453     // TOBTEC: look for rec-seeds with hits in tob1 and 2 only.
0454     int tot_cnt = 0, no_mkf_cnt = 0;
0455 
0456     for (auto& [l, simt_ptr]: m_sim_map)
0457     {
0458       auto &simt = * simt_ptr;
0459       auto &ckft = * m_ckf_map[l];
0460       auto mi = m_mkf_map.find(l);
0461 
0462       bool print_all = print_all_def;
0463 
0464       // TOBTEC: look for rec-seeds with hits in tob1 and 2 only.
0465       bool select = true;
0466       {
0467         auto &ckf_seed = ev.seedTracks_[ckft.label()];
0468         for (int hi = 0; hi < ckf_seed.nTotalHits(); ++hi) {
0469           const HitOnTrack hot = ckf_seed.getHitOnTrack(hi);
0470           if (hot.index >= 0 && (hot.layer < 10 || hot.layer > 13)) {
0471             select = false;
0472             break;
0473           }
0474         }
0475       }
0476       if ( ! select) continue;
0477 
0478       ++tot_cnt;
0479       //print_all = true;
0480 
0481       if (mi != m_mkf_map.end())
0482       {
0483         auto &mkft = * mi->second;
0484         mkf_cnt++;
0485         if (mkft.nFoundHits() < ckft.nFoundHits()) ++less_hits;
0486         if (mkft.nFoundHits() > ckft.nFoundHits()) ++more_hits;
0487 
0488         CheckMkFitLayerPlanVsReferenceHits(mkft, ckft, "CKF");
0489         // CheckMkFitLayerPlanVsReferenceHits(mkft, simt, "SIM");
0490 
0491         (void) print_all;
0492         if (/* itconf.m_track_algorithm == 10 ||*/ print_all) {
0493           // ckf label is wrong when validation is on (even quality val) for mixedTriplet, pixelless and tobtec
0494           // as seed tracks get removed for non-mkfit iterations and indices from rec-tracks are no longer valid.
0495           auto &ckf_seed = ev.seedTracks_[ckft.label()];
0496           auto &mkf_seed = m_seeds[mkft.label()];
0497           print("ckf  ", 0, ckft, ev);
0498           print("mkfit", 0, mkft, ev);
0499           print("sim  ", 0, simt, ev);
0500           print("ckf seed", 0, ckf_seed, ev);
0501           print("mkf seed", 0, mkf_seed, ev);
0502           printf("------------------------------------------------------\n");
0503 
0504           TrackVec ssss;
0505           ssss.push_back(mkf_seed);
0506 
0507           IterationSeedPartition pppp(1);
0508           IterationConfig::get_seed_partitioner("phase1:1:debug")(Config::TrkInfo, ssss, *m_eoh, pppp);
0509 
0510           printf("------------------------------------------------------\n");
0511           printf("\n");
0512         }
0513       }
0514       else
0515       {
0516         printf("\n!!!!! No mkfit track with this label.\n\n");
0517         ++no_mkf_cnt;
0518 
0519         auto &ckf_seed = ev.seedTracks_[ckft.label()];
0520         print("ckf ", 0, ckft, ev);
0521         print("sim ", 0, simt, ev);
0522         print("ckf seed", 0, ckf_seed, ev);
0523         auto smi = m_seed_map.find(l);
0524         if (smi != m_seed_map.end())
0525           print("seed with matching label", 0, *smi->second, ev);
0526         printf("------------------------------------------------------\n");
0527       }
0528     }
0529 
0530     printf("mkFit found %d, matching_label=%d, less_hits=%d, more_hits=%d  [algo=%d (%s)]\n",
0531            (int) ev.fitTracks_.size(), mkf_cnt, less_hits, more_hits,
0532            itconf.m_track_algorithm, TrackBase::algoint_to_cstr(itconf.m_track_algorithm));
0533 
0534     if (tot_cnt > 0) {
0535       printf("\ntobtec tob1/2 tot=%d no_mkf=%d (%f%%)\n",
0536             tot_cnt, no_mkf_cnt, 100.0 * no_mkf_cnt / tot_cnt);
0537     } else {
0538       printf("\nNo CKF tracks with seed hits in TOB1/2 found (need iteration idx 8, TobTec?)\n");
0539     }
0540 
0541     printf("-------------------------------------------------------------------------------------------\n");
0542     printf("-------------------------------------------------------------------------------------------\n");
0543     printf("\n");
0544   }
0545 
0546 }