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     }
0131
0132     void Quality::quality_reset() { m_cnt = m_cnt1 = m_cnt2 = m_cnt_8 = m_cnt1_8 = m_cnt2_8 = m_cnt_nomc = 0; }
0133
0134     void Quality::quality_process(Event *event, Track &tkcand, const int itrack, std::map<int, int> &cmsswLabelToPos) {
0135       // KPM: Do not use this method for validating CMSSW tracks if we ever build a DumbCMSSW function for them to print out...
0136       // as we would need to access seeds through map of seed ids...
0137
0138       // initialize track extra (input original seed label)
0139       const auto label = tkcand.label();
0140       TrackExtra extra(label);
0141
0142       // track_print(event, tkcand, "quality_process -> track_print:");
0143
0144       // access temp seed trk and set matching seed hits
0145       const auto &seed = event->seedTracks_[itrack];
0146       extra.findMatchingSeedHits(tkcand, seed, event->layerHits_);
0147
0148       // set mcTrackID through 50% hit matching after seed
0149       extra.setMCTrackIDInfo(
0150           tkcand, event->layerHits_, event->simHitsInfo_, event->simTracks_, false, (Config::seedInput == simSeeds));
0151       const int mctrk = extra.mcTrackID();
0152
0153       //  int mctrk = tkcand.label(); // assumes 100% "efficiency"
0154
0155       const float pT = tkcand.pT();
0156       float pTmc = 0.f, etamc = 0.f, phimc = 0.f;
0157       float pTr = 0.f;
0158       int nfoundmc = -1;
0159
0160       if (mctrk < 0 || static_cast<size_t>(mctrk) >= event->simTracks_.size()) {
0161         ++m_cnt_nomc;
0162         dprint("XX bad track idx " << mctrk << ", orig label was " << label);
0163       } else {
0164         auto &simtrack = event->simTracks_[mctrk];
0165         pTmc = simtrack.pT();
0166         etamc = simtrack.momEta();
0167         phimc = simtrack.momPhi();
0168         pTr = pT / pTmc;
0169
0170         nfoundmc = simtrack.nUniqueLayers();
0171
0172         ++m_cnt;
0173         if (pTr > 0.9 && pTr < 1.1)
0174           ++m_cnt1;
0175         if (pTr > 0.8 && pTr < 1.2)
0176           ++m_cnt2;
0177
0178         if (tkcand.nFoundHits() >= 0.8f * nfoundmc) {
0179           ++m_cnt_8;
0180           if (pTr > 0.9 && pTr < 1.1)
0181             ++m_cnt1_8;
0182           if (pTr > 0.8 && pTr < 1.2)
0183             ++m_cnt2_8;
0184         }
0185
0186         // perl -ne 'print if m/FOUND_LABEL\s+[-\d]+/o;' | sort -k2 -n
0187         // grep "FOUND_LABEL" | sort -n -k 8,8 -k 2,2
0188         // printf("FOUND_LABEL %6d  pT_mc= %8.2f eta_mc= %8.2f event= %d\n", label, pTmc, etamc, event->evtID());
0189       }
0190
0191       float pTcmssw = 0.f, etacmssw = 0.f, phicmssw = 0.f;
0192       int nfoundcmssw = -1;
0193       if (Config::dumpForPlots && Config::readCmsswTracks) {
0194         if (cmsswLabelToPos.count(label)) {
0195           auto &cmsswtrack = event->cmsswTracks_[cmsswLabelToPos[label]];
0196           pTcmssw = cmsswtrack.pT();
0197           etacmssw = cmsswtrack.momEta();
0198           phicmssw = cmsswtrack.swimPhiToR(tkcand.x(), tkcand.y());  // to get rough estimate of diff in phi
0199           nfoundcmssw = cmsswtrack.nUniqueLayers();
0200         }
0201       }
0202
0203       if (!Config::silent && Config::dumpForPlots) {
0204         std::lock_guard<std::mutex> printlock(Event::printmutex);
0205         printf(
0206             "MX - found track with chi2= %6.3f nFoundHits= %2d pT= %7.4f eta= %7.4f phi= %7.4f nfoundmc= %2d pTmc= "
0207             "%7.4f etamc= %7.4f phimc= %7.4f nfoundcmssw= %2d pTcmssw= %7.4f etacmssw= %7.4f phicmssw= %7.4f lab= %d\n",
0208             tkcand.chi2(),
0209             tkcand.nFoundHits(),
0210             pT,
0211             tkcand.momEta(),
0212             tkcand.momPhi(),
0213             nfoundmc,
0214             pTmc,
0215             etamc,
0216             phimc,
0217             nfoundcmssw,
0218             pTcmssw,
0219             etacmssw,
0220             phicmssw,
0221             label);
0222       }
0223     }
0224
0225     void Quality::quality_print() {
0226       if (!Config::silent) {
0227         std::lock_guard<std::mutex> printlock(Event::printmutex);
0228         std::cout << "found tracks=" << m_cnt << "  in pT 10%=" << m_cnt1 << "  in pT 20%=" << m_cnt2
0229                   << "     no_mc_assoc=" << m_cnt_nomc << std::endl;
0230         std::cout << "  nH >= 80% =" << m_cnt_8 << "  in pT 10%=" << m_cnt1_8 << "  in pT 20%=" << m_cnt2_8
0231                   << std::endl;
0232       }
0233     }
0234
0235     //------------------------------------------------------------------------------
0236     // Root validation
0237     //------------------------------------------------------------------------------
0238
0239     void root_val_dumb_cmssw(Event *event) {
0240       // get labels correct first
0241       event->relabel_bad_seedtracks();
0242       event->relabel_cmsswtracks_from_seeds();
0243
0244       //collection cleaning
0245       if (Config::nItersCMSSW > 0)
0246         event->select_tracks_iter(Config::nItersCMSSW);
0247
0248       // set the track collections to each other
0249       event->candidateTracks_ = event->cmsswTracks_;
0250       event->fitTracks_ = event->candidateTracks_;
0251
0252       // prep the tracks + extras
0253       prep_simtracks(event);
0254       prep_recotracks(event);
0255
0256       // validate
0257       event->validate();
0258     }
0259
0260     void root_val(Event *event) {
0261       // score the tracks
0262       score_tracks(event->seedTracks_);
0263       score_tracks(event->candidateTracks_);
0264
0265       // deal with fit tracks
0266       if (Config::backwardFit) {
0267         score_tracks(event->fitTracks_);
0268       } else
0269         event->fitTracks_ = event->candidateTracks_;
0270
0271       // sort hits + make extras, align if needed
0272       prep_recotracks(event);
0273       if (Config::cmssw_val)
0274         prep_cmsswtracks(event);
0275
0276       // validate
0277       event->validate();
0278     }
0279
0280     void prep_recotracks(Event *event) {
0281       // seed tracks extras always needed
0282       if (Config::sim_val || Config::sim_val_for_cmssw) {
0283         prep_tracks(event, event->seedTracks_, event->seedTracksExtra_, true);
0284       } else if (Config::cmssw_val)  // seed tracks are not validated, labels used for maps --> do NOT align index and labels!
0285       {
0286         prep_tracks(event, event->seedTracks_, event->seedTracksExtra_, false);
0287       }
0288
0289       // make extras + align index == label() for candidate tracks
0290       prep_tracks(event, event->candidateTracks_, event->candidateTracksExtra_, true);
0291       prep_tracks(event, event->fitTracks_, event->fitTracksExtra_, true);
0292     }
0293
0294     void prep_simtracks(Event *event) {
0295       // First prep sim tracks to have hits sorted, then mark unfindable if too short
0296       prep_reftracks(event, event->simTracks_, event->simTracksExtra_, false);
0297
0298       // Now, make sure sim track shares at least four hits with a single cmssw seed.
0299       // This ensures we factor out any weakness from CMSSW
0300
0301       // First, make a make a map of [lyr][hit idx].vector(seed trk labels)
0302       LayIdxIDVecMapMap seedHitIDMap;
0303       std::map<int, int> labelNHitsMap;
0304       std::map<int, int> labelAlgoMap;
0305       std::map<int, std::vector<int>> labelSeedHitsMap;
0306       for (const auto &seedtrack : event->seedTracks_) {
0307         for (int ihit = 0; ihit < seedtrack.nTotalHits(); ihit++) {
0308           const auto lyr = seedtrack.getHitLyr(ihit);
0309           const auto idx = seedtrack.getHitIdx(ihit);
0310
0311           if (lyr < 0 || idx < 0)
0312             continue;  // standard check
0313           seedHitIDMap[lyr][idx].push_back(seedtrack.label());
0314           labelSeedHitsMap[seedtrack.label()].push_back(lyr);
0315         }
0316         labelNHitsMap[seedtrack.label()] = seedtrack.nTotalHits();
0317         labelAlgoMap[seedtrack.label()] = seedtrack.algoint();
0318       }
0319
0320       // Then, loop over sim tracks, and add up how many lyrs they possess of a single seed track
0321       unsigned int count = 0;
0322       for (auto &simtrack : event->simTracks_) {
0323         if (simtrack.isNotFindable())
0324           continue;  // skip ones we already know are bad
0325         TrkIDLaySetMap seedIDMap;
0326         for (int ihit = 0; ihit < simtrack.nTotalHits(); ihit++) {
0327           const auto lyr = simtrack.getHitLyr(ihit);
0328           const auto idx = simtrack.getHitIdx(ihit);
0329
0330           if (lyr < 0 || idx < 0)
0331             continue;  // standard check
0332
0333           if (!seedHitIDMap.count(lyr))
0334             continue;  // ensure seed hit map has at least one entry for this layer
0335           if (!seedHitIDMap.at(lyr).count(idx))
0336             continue;  // ensure seed hit map has at least one entry for this idx
0337
0338           for (const auto label : seedHitIDMap.at(lyr).at(idx)) {
0339             const auto &seedLayers = labelSeedHitsMap[label];
0340             if (std::find(seedLayers.begin(), seedLayers.end(), lyr) != seedLayers.end())  //seed check moved here
0341               seedIDMap[label].emplace(lyr);
0342           }
0343         }
0344
0345         // now see if one of the seedIDs matched has at least 4 hits!
0346         bool isSimSeed = false;
0347         for (const auto &seedIDpair : seedIDMap) {
0348           if ((int)seedIDpair.second.size() == labelNHitsMap[seedIDpair.first]) {
0349             isSimSeed = true;
0350             if (Config::mtvRequireSeeds)
0351               simtrack.setAlgoint(labelAlgoMap[seedIDpair.first]);
0352             if (Config::mtvRequireSeeds)
0353               event->simTracksExtra_[count].addAlgo(labelAlgoMap[seedIDpair.first]);
0354             //break;
0355           }
0356         }
0357         if (Config::mtvLikeValidation) {
0358           // Apply MTV selection criteria and then return
0359           if (simtrack.prodType() != Track::ProdType::Signal || simtrack.charge() == 0 || simtrack.posR() > 2.5 ||
0360               std::abs(simtrack.z()) > 30 || std::abs(simtrack.momEta()) > 3.0)
0361             simtrack.setNotFindable();
0362           else if (Config::mtvRequireSeeds && !isSimSeed)
0363             simtrack.setNotFindable();
0364         } else {
0365           // set findability based on bool isSimSeed
0366           if (!isSimSeed)
0367             simtrack.setNotFindable();
0368         }
0369         count++;
0370       }
0371     }
0372
0373     void prep_cmsswtracks(Event *event) { prep_reftracks(event, event->cmsswTracks_, event->cmsswTracksExtra_, true); }
0374
0375     void prep_reftracks(Event *event, TrackVec &tracks, TrackExtraVec &extras, const bool realigntracks) {
0376       prep_tracks(event, tracks, extras, realigntracks);
0377
0378       // mark cmsswtracks as unfindable if too short
0379       for (auto &track : tracks) {
0380         const int nlyr = track.nUniqueLayers();
0381         if (nlyr < Config::cmsSelMinLayers)
0382           track.setNotFindable();
0383       }
0384     }
0385
0386     void prep_tracks(Event *event, TrackVec &tracks, TrackExtraVec &extras, const bool realigntracks) {
0387       for (size_t i = 0; i < tracks.size(); i++) {
0388         extras.emplace_back(tracks[i].label());
0389       }
0390       if (realigntracks)
0391         event->validation_.alignTracks(tracks, extras, false);
0392     }
0393
0394     void score_tracks(TrackVec &tracks) {
0395       auto score_func = IterationConfig::get_track_scorer("default");
0396       for (auto &track : tracks) {
0397         track.setScore(getScoreCand(score_func, track));
0398       }
0399     }
0400
0401   }  // namespace StdSeq
0402
0403 }  // namespace mkfit