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