MkFitCMS/standalone/MkStandaloneSeqs.cc

0001 #include "RecoTracker/MkFitCMS/standalone/MkStandaloneSeqs.h"
0002 #include "RecoTracker/MkFitCMS/interface/MkStdSeqs.h"
0003
0004 #include "RecoTracker/MkFitCore/interface/HitStructures.h"
0005 #include "RecoTracker/MkFitCore/standalone/Event.h"
0006 #include "RecoTracker/MkFitCore/interface/IterationConfig.h"
0007
0008 #include "RecoTracker/MkFitCore/src/Debug.h"
0009
0010 #include "oneapi/tbb/parallel_for.h"
0011
0012 namespace mkfit {
0013
0014   namespace StdSeq {
0015
0016     //=========================================================================
0017     // Hit processing
0018     //=========================================================================
0019
0020     void loadHitsAndBeamSpot(Event &ev, EventOfHits &eoh) {
0021       eoh.reset();
0022
0023       // fill vector of hits in each layer
0024       // XXXXMT: Does it really makes sense to multi-thread this?
0025       TBB_PARALLEL_FOR(tbb::blocked_range<int>(0, ev.layerHits_.size()), [&](const tbb::blocked_range<int> &layers) {
0026         for (int ilay = layers.begin(); ilay < layers.end(); ++ilay) {
0027           eoh.suckInHits(ilay, ev.layerHits_[ilay]);
0028         }
0029       });
0030       eoh.setBeamSpot(ev.beamSpot_);
0031     }
0032
0033     void handle_duplicates(Event *) {
0034       /*
0035       // Mark tracks as duplicates; if within CMSSW, remove duplicate tracks from fit or candidate track collection
0036       if (Config::removeDuplicates) {
0037         if (Config::quality_val || Config::sim_val || Config::cmssw_val) {
0038           clean_duplicates(event->candidateTracks_);
0039           if (Config::backwardFit)
0040             clean_duplicates(event->fitTracks_);
0041         }
0042         // For the MEIF benchmarks and the stress tests, no validation flags are set so we will enter this block
0043         else {
0044           // Only care about the candidate tracks here; no need to run the duplicate removal on both candidate and fit tracks
0045           clean_duplicates(event->candidateTracks_);
0046         }
0047       }
0048       */
0049     }
0050
0051     //=========================================================================
0052     // Random stuff
0053     //=========================================================================
0054
0055     void dump_simtracks(Event *event) {
0056       // Ripped out of MkBuilder::begin_event, ifdefed under DEBUG
0057
0058       std::vector<Track> &simtracks = event->simTracks_;
0059
0060       for (int itrack = 0; itrack < (int)simtracks.size(); ++itrack) {
0061         // bool debug = true;
0062         Track track = simtracks[itrack];
0063         // simtracks are initially written with label = index; uncomment in case tracks were edited
0064         // if (track.label() != itrack) {
0065         //   dprintf("Bad label for simtrack %d -- %d\n", itrack, track.label());
0066         // }
0067
0068         dprint("MX - simtrack with nHits=" << track.nFoundHits() << " chi2=" << track.chi2() << " pT=" << track.pT()
0069                                            << " phi=" << track.momPhi() << " eta=" << track.momEta());
0070       }
0071
0072       for (int itrack = 0; itrack < (int)simtracks.size(); ++itrack) {
0073         for (int ihit = 0; ihit < simtracks[itrack].nFoundHits(); ++ihit) {
0074           dprint("track #" << itrack << " hit #" << ihit
0075                            << " hit pos=" << simtracks[itrack].hitsVector(event->layerHits_)[ihit].position()
0076                            << " phi=" << simtracks[itrack].hitsVector(event->layerHits_)[ihit].phi());
0077         }
0078       }
0079     }
0080
0081     void track_print(Event *event, const Track &t, const char *pref) {
0082       printf("%s with q=%+i pT=%7.3f eta=% 7.3f nHits=%2d  label=%4d\nState:\n",
0083              pref,
0084              t.charge(),
0085              t.pT(),
0086              t.momEta(),
0087              t.nFoundHits(),
0088              t.label());
0089
0090       print(t.state());
0091
0092       printf("Hits:\n");
0093       for (int ih = 0; ih < t.nTotalHits(); ++ih) {
0094         int lyr = t.getHitLyr(ih);
0095         int idx = t.getHitIdx(ih);
0096         if (idx >= 0) {
0097           const Hit &hit = event->layerHits_[lyr][idx];
0098           printf("    hit %2d lyr=%2d idx=%4d pos r=%7.3f z=% 8.3f   mc_hit=%4d mc_trk=%4d\n",
0099                  ih,
0100                  lyr,
0101                  idx,
0102                  hit.r(),
0103                  hit.z(),
0104                  hit.mcHitID(),
0105                  hit.mcTrackID(event->simHitsInfo_));
0106         } else
0107           printf("    hit %2d        idx=%i\n", ih, t.getHitIdx(ih));
0108       }
0109     }
0110
0111     //------------------------------------------------------------------------------
0112     // Non-ROOT validation
0113     //------------------------------------------------------------------------------
0114
0115     void Quality::quality_val(Event *event) {
0116       quality_reset();
0117
0118       std::map<int, int> cmsswLabelToPos;
0119       if (Config::dumpForPlots && Config::readCmsswTracks) {
0120         for (size_t itrack = 0; itrack < event->cmsswTracks_.size(); itrack++) {
0121           cmsswLabelToPos[event->cmsswTracks_[itrack].label()] = itrack;
0122         }
0123       }
0124
0125       for (size_t itrack = 0; itrack < event->candidateTracks_.size(); itrack++) {
0126         quality_process(event, event->candidateTracks_[itrack], itrack, cmsswLabelToPos);
0127       }
0128
0129       quality_print();
0130       add_to_quality_sum(*this);
0131     }
0132
0133     void Quality::quality_reset() { m_cnt = m_cnt1 = m_cnt2 = m_cnt_8 = m_cnt1_8 = m_cnt2_8 = m_cnt_nomc = 0; }
0134
0135     void Quality::quality_process(Event *event, Track &tkcand, const int itrack, std::map<int, int> &cmsswLabelToPos) {
0136       // KPM: Do not use this method for validating CMSSW tracks if we ever build a DumbCMSSW function for them to print out...
0137       // as we would need to access seeds through map of seed ids...
0138
0139       // initialize track extra (input original seed label)
0140       const auto label = tkcand.label();
0141       TrackExtra extra(label);
0142
0143       // track_print(event, tkcand, "quality_process -> track_print:");
0144
0145       // access temp seed trk and set matching seed hits
0146       const auto &seed = event->seedTracks_[itrack];
0147       extra.findMatchingSeedHits(tkcand, seed, event->layerHits_);
0148
0149       // set mcTrackID through 50% hit matching after seed
0150       extra.setMCTrackIDInfo(
0151           tkcand, event->layerHits_, event->simHitsInfo_, event->simTracks_, false, (Config::seedInput == simSeeds));
0152       const int mctrk = extra.mcTrackID();
0153
0154       //  int mctrk = tkcand.label(); // assumes 100% "efficiency"
0155
0156       const float pT = tkcand.pT();
0157       float pTmc = 0.f, etamc = 0.f, phimc = 0.f;
0158       float pTr = 0.f;
0159       int nfoundmc = -1;
0160
0161       if (mctrk < 0 || static_cast<size_t>(mctrk) >= event->simTracks_.size()) {
0162         ++m_cnt_nomc;
0163         dprintf("XX bad track idx %d, orig label was %d\n", mctrk, label);
0164       } else {
0165         auto &simtrack = event->simTracks_[mctrk];
0166         pTmc = simtrack.pT();
0167         etamc = simtrack.momEta();
0168         phimc = simtrack.momPhi();
0169         pTr = pT / pTmc;
0170
0171         nfoundmc = simtrack.nUniqueLayers();
0172
0173         ++m_cnt;
0174         if (pTr > 0.9 && pTr < 1.1)
0175           ++m_cnt1;
0176         if (pTr > 0.8 && pTr < 1.2)
0177           ++m_cnt2;
0178
0179         if (tkcand.nFoundHits() >= 0.8f * nfoundmc) {
0180           ++m_cnt_8;
0181           if (pTr > 0.9 && pTr < 1.1)
0182             ++m_cnt1_8;
0183           if (pTr > 0.8 && pTr < 1.2)
0184             ++m_cnt2_8;
0185         }
0186
0187         // perl -ne 'print if m/FOUND_LABEL\s+[-\d]+/o;' | sort -k2 -n
0188         // grep "FOUND_LABEL" | sort -n -k 8,8 -k 2,2
0189         // printf("FOUND_LABEL %6d  pT_mc= %8.2f eta_mc= %8.2f event= %d\n", label, pTmc, etamc, event->evtID());
0190       }
0191
0192       float pTcmssw = 0.f, etacmssw = 0.f, phicmssw = 0.f;
0193       int nfoundcmssw = -1;
0194       if (Config::dumpForPlots && Config::readCmsswTracks) {
0195         if (cmsswLabelToPos.count(label)) {
0196           auto &cmsswtrack = event->cmsswTracks_[cmsswLabelToPos[label]];
0197           pTcmssw = cmsswtrack.pT();
0198           etacmssw = cmsswtrack.momEta();
0199           phicmssw = cmsswtrack.swimPhiToR(tkcand.x(), tkcand.y());  // to get rough estimate of diff in phi
0200           nfoundcmssw = cmsswtrack.nUniqueLayers();
0201         }
0202       }
0203
0204       if (!Config::silent && Config::dumpForPlots) {
0205         std::lock_guard<std::mutex> printlock(Event::printmutex);
0206         printf(
0207             "MX - found track with chi2= %6.3f nFoundHits= %2d pT= %7.4f eta= %7.4f phi= %7.4f nfoundmc= %2d pTmc= "
0208             "%7.4f etamc= %7.4f phimc= %7.4f nfoundcmssw= %2d pTcmssw= %7.4f etacmssw= %7.4f phicmssw= %7.4f lab= %d\n",
0209             tkcand.chi2(),
0210             tkcand.nFoundHits(),
0211             pT,
0212             tkcand.momEta(),
0213             tkcand.momPhi(),
0214             nfoundmc,
0215             pTmc,
0216             etamc,
0217             phimc,
0218             nfoundcmssw,
0219             pTcmssw,
0220             etacmssw,
0221             phicmssw,
0222             label);
0223       }
0224     }
0225
0226     void Quality::quality_print() {
0227       if (!Config::silent) {
0228         std::lock_guard<std::mutex> printlock(Event::printmutex);
0229         std::cout << "found tracks=" << m_cnt << "  in pT 10%=" << m_cnt1 << "  in pT 20%=" << m_cnt2
0230                   << "     no_mc_assoc=" << m_cnt_nomc << std::endl;
0231         std::cout << "  nH >= 80% =" << m_cnt_8 << "  in pT 10%=" << m_cnt1_8 << "  in pT 20%=" << m_cnt2_8
0232                   << std::endl;
0233       }
0234     }
0235
0236     Quality Quality::s_quality_sum;
0237
0238     void Quality::add_to_quality_sum(const Quality &q) {
0239       static std::mutex q_mutex;
0240       std::lock_guard<std::mutex> q_lock(q_mutex);
0241
0242       s_quality_sum.m_cnt += q.m_cnt;
0243       s_quality_sum.m_cnt1 += q.m_cnt1;
0244       s_quality_sum.m_cnt2 += q.m_cnt2;
0245       s_quality_sum.m_cnt_8 += q.m_cnt_8;
0246       s_quality_sum.m_cnt1_8 += q.m_cnt1_8;
0247       s_quality_sum.m_cnt2_8 += q.m_cnt2_8;
0248       s_quality_sum.m_cnt_nomc += q.m_cnt_nomc;
0249     }
0250
0251     //------------------------------------------------------------------------------
0252     // Root validation
0253     //------------------------------------------------------------------------------
0254
0255     void root_val_dumb_cmssw(Event *event) {
0256       // get labels correct first
0257       event->relabel_bad_seedtracks();
0258       event->relabel_cmsswtracks_from_seeds();
0259
0260       //collection cleaning
0261       if (Config::nItersCMSSW > 0)
0262         event->select_tracks_iter(Config::nItersCMSSW);
0263
0264       // set the track collections to each other
0265       event->candidateTracks_ = event->cmsswTracks_;
0266       event->fitTracks_ = event->candidateTracks_;
0267
0268       // prep the tracks + extras
0269       prep_simtracks(event);
0270       prep_recotracks(event);
0271
0272       // validate
0273       event->validate();
0274     }
0275
0276     void root_val(Event *event) {
0277       // score the tracks
0278       score_tracks(event->seedTracks_);
0279       score_tracks(event->candidateTracks_);
0280
0281       // deal with fit tracks
0282       if (Config::backwardFit) {
0283         score_tracks(event->fitTracks_);
0284       } else
0285         event->fitTracks_ = event->candidateTracks_;
0286
0287       // sort hits + make extras, align if needed
0288       prep_recotracks(event);
0289       if (Config::cmssw_val)
0290         prep_cmsswtracks(event);
0291
0292       // validate
0293       event->validate();
0294     }
0295
0296     void prep_recotracks(Event *event) {
0297       // seed tracks extras always needed
0298       if (Config::sim_val || Config::sim_val_for_cmssw) {
0299         prep_tracks(event, event->seedTracks_, event->seedTracksExtra_, true);
0300       } else if (Config::cmssw_val)  // seed tracks are not validated, labels used for maps --> do NOT align index and labels!
0301       {
0302         prep_tracks(event, event->seedTracks_, event->seedTracksExtra_, false);
0303       }
0304
0305       // make extras + align index == label() for candidate tracks
0306       prep_tracks(event, event->candidateTracks_, event->candidateTracksExtra_, true);
0307       prep_tracks(event, event->fitTracks_, event->fitTracksExtra_, true);
0308     }
0309
0310     void prep_simtracks(Event *event) {
0311       // First prep sim tracks to have hits sorted, then mark unfindable if too short
0312       prep_reftracks(event, event->simTracks_, event->simTracksExtra_, false);
0313
0314       // Now, make sure sim track shares at least four hits with a single cmssw seed.
0315       // This ensures we factor out any weakness from CMSSW
0316
0317       // First, make a make a map of [lyr][hit idx].vector(seed trk labels)
0318       LayIdxIDVecMapMap seedHitIDMap;
0319       std::map<int, int> labelNHitsMap;
0320       std::map<int, int> labelAlgoMap;
0321       std::map<int, std::vector<int>> labelSeedHitsMap;
0322       for (const auto &seedtrack : event->seedTracks_) {
0323         for (int ihit = 0; ihit < seedtrack.nTotalHits(); ihit++) {
0324           const auto lyr = seedtrack.getHitLyr(ihit);
0325           const auto idx = seedtrack.getHitIdx(ihit);
0326
0327           if (lyr < 0 || idx < 0)
0328             continue;  // standard check
0329           seedHitIDMap[lyr][idx].push_back(seedtrack.label());
0330           labelSeedHitsMap[seedtrack.label()].push_back(lyr);
0331         }
0332         labelNHitsMap[seedtrack.label()] = seedtrack.nTotalHits();
0333         labelAlgoMap[seedtrack.label()] = seedtrack.algoint();
0334       }
0335
0336       // Then, loop over sim tracks, and add up how many lyrs they possess of a single seed track
0337       unsigned int count = 0;
0338       for (auto &simtrack : event->simTracks_) {
0339         if (simtrack.isNotFindable())
0340           continue;  // skip ones we already know are bad
0341         TrkIDLaySetMap seedIDMap;
0342         for (int ihit = 0; ihit < simtrack.nTotalHits(); ihit++) {
0343           const auto lyr = simtrack.getHitLyr(ihit);
0344           const auto idx = simtrack.getHitIdx(ihit);
0345
0346           if (lyr < 0 || idx < 0)
0347             continue;  // standard check
0348
0349           if (!seedHitIDMap.count(lyr))
0350             continue;  // ensure seed hit map has at least one entry for this layer
0351           if (!seedHitIDMap.at(lyr).count(idx))
0352             continue;  // ensure seed hit map has at least one entry for this idx
0353
0354           for (const auto label : seedHitIDMap.at(lyr).at(idx)) {
0355             const auto &seedLayers = labelSeedHitsMap[label];
0356             if (std::find(seedLayers.begin(), seedLayers.end(), lyr) != seedLayers.end())  //seed check moved here
0357               seedIDMap[label].emplace(lyr);
0358           }
0359         }
0360
0361         // now see if one of the seedIDs matched has at least 4 hits!
0362         bool isSimSeed = false;
0363         for (const auto &seedIDpair : seedIDMap) {
0364           if ((int)seedIDpair.second.size() == labelNHitsMap[seedIDpair.first]) {
0365             isSimSeed = true;
0366             if (Config::mtvRequireSeeds)
0367               simtrack.setAlgoint(labelAlgoMap[seedIDpair.first]);
0368             if (Config::mtvRequireSeeds)
0369               event->simTracksExtra_[count].addAlgo(labelAlgoMap[seedIDpair.first]);
0370             //break;
0371           }
0372         }
0373         if (Config::mtvLikeValidation) {
0374           // Apply MTV selection criteria and then return
0375           if (simtrack.prodType() != Track::ProdType::Signal || simtrack.charge() == 0 || simtrack.posR() > 2.5 ||
0376               std::abs(simtrack.z()) > 30 || std::abs(simtrack.momEta()) > 3.0)
0377             simtrack.setNotFindable();
0378           else if (Config::mtvRequireSeeds && !isSimSeed)
0379             simtrack.setNotFindable();
0380         } else {
0381           // set findability based on bool isSimSeed
0382           if (!isSimSeed)
0383             simtrack.setNotFindable();
0384         }
0385         count++;
0386       }
0387     }
0388
0389     void prep_cmsswtracks(Event *event) { prep_reftracks(event, event->cmsswTracks_, event->cmsswTracksExtra_, true); }
0390
0391     void prep_reftracks(Event *event, TrackVec &tracks, TrackExtraVec &extras, const bool realigntracks) {
0392       prep_tracks(event, tracks, extras, realigntracks);
0393
0394       // mark cmsswtracks as unfindable if too short
0395       for (auto &track : tracks) {
0396         const int nlyr = track.nUniqueLayers();
0397         if (nlyr < Config::cmsSelMinLayers)
0398           track.setNotFindable();
0399       }
0400     }
0401
0402     void prep_tracks(Event *event, TrackVec &tracks, TrackExtraVec &extras, const bool realigntracks) {
0403       for (size_t i = 0; i < tracks.size(); i++) {
0404         extras.emplace_back(tracks[i].label());
0405       }
0406       if (realigntracks)
0407         event->validation_.alignTracks(tracks, extras, false);
0408     }
0409
0410     void score_tracks(TrackVec &tracks) {
0411       auto score_func = IterationConfig::get_track_scorer("default");
0412       for (auto &track : tracks) {
0413         track.setScore(getScoreCand(score_func, track));
0414       }
0415     }
0416
0417   }  // namespace StdSeq
0418
0419 }  // namespace mkfit