File indexing completed on 2022-02-21 23:14:11
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
0017
0018
0019 void loadHitsAndBeamSpot(Event &ev, EventOfHits &eoh) {
0020 eoh.reset();
0021
0022
0023
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
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
0041 else {
0042
0043 find_duplicates(event->candidateTracks_);
0044 }
0045 }
0046 }
0047
0048
0049
0050
0051
0052 void dump_simtracks(Event *event) {
0053
0054
0055 std::vector<Track> &simtracks = event->simTracks_;
0056
0057 for (int itrack = 0; itrack < (int)simtracks.size(); ++itrack) {
0058
0059 Track track = simtracks[itrack];
0060
0061
0062
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
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
0133
0134
0135
0136 const auto label = tkcand.label();
0137 TrackExtra extra(label);
0138
0139
0140
0141
0142 const auto &seed = event->seedTracks_[itrack];
0143 extra.findMatchingSeedHits(tkcand, seed, event->layerHits_);
0144
0145
0146 extra.setMCTrackIDInfo(
0147 tkcand, event->layerHits_, event->simHitsInfo_, event->simTracks_, false, (Config::seedInput == simSeeds));
0148 const int mctrk = extra.mcTrackID();
0149
0150
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
0184
0185
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());
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
0239
0240
0241 void root_val_dumb_cmssw(Event *event) {
0242
0243 event->relabel_bad_seedtracks();
0244 event->relabel_cmsswtracks_from_seeds();
0245
0246
0247 if (Config::nItersCMSSW > 0)
0248 event->select_tracks_iter(Config::nItersCMSSW);
0249
0250
0251 event->candidateTracks_ = event->cmsswTracks_;
0252 event->fitTracks_ = event->candidateTracks_;
0253
0254
0255 prep_simtracks(event);
0256 prep_recotracks(event);
0257
0258
0259 event->validate();
0260 }
0261
0262 void root_val(Event *event) {
0263
0264 score_tracks(event->seedTracks_);
0265 score_tracks(event->candidateTracks_);
0266
0267
0268 if (Config::backwardFit) {
0269 score_tracks(event->fitTracks_);
0270 } else
0271 event->fitTracks_ = event->candidateTracks_;
0272
0273
0274 prep_recotracks(event);
0275 if (Config::cmssw_val)
0276 prep_cmsswtracks(event);
0277
0278
0279 event->validate();
0280 }
0281
0282 void prep_recotracks(Event *event) {
0283
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)
0287 {
0288 prep_tracks(event, event->seedTracks_, event->seedTracksExtra_, false);
0289 }
0290
0291
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
0298 prep_reftracks(event, event->simTracks_, event->simTracksExtra_, false);
0299
0300
0301
0302
0303
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;
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
0323 unsigned int count = 0;
0324 for (auto &simtrack : event->simTracks_) {
0325 if (simtrack.isNotFindable())
0326 continue;
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;
0334
0335 if (!seedHitIDMap.count(lyr))
0336 continue;
0337 if (!seedHitIDMap.at(lyr).count(idx))
0338 continue;
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())
0343 seedIDMap[label].emplace(lyr);
0344 }
0345 }
0346
0347
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
0357 }
0358 }
0359 if (Config::mtvLikeValidation) {
0360
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
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
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 }
0403
0404 }