File indexing completed on 2024-04-06 12:28:25
0001 #include "RecoTracker/MkFitCore/standalone/TrackExtra.h"
0002 #include "RecoTracker/MkFitCore/standalone/ConfigStandalone.h"
0003
0004
0005 #include "RecoTracker/MkFitCore/src/Debug.h"
0006
0007 namespace mkfit {
0008
0009
0010
0011
0012
0013 void TrackExtra::findMatchingSeedHits(const Track& reco_trk,
0014 const Track& seed_trk,
0015 const std::vector<HitVec>& layerHits) {
0016
0017 for (int reco_ihit = 0; reco_ihit < reco_trk.nTotalHits(); ++reco_ihit) {
0018 const int reco_lyr = reco_trk.getHitLyr(reco_ihit);
0019 const int reco_idx = reco_trk.getHitIdx(reco_ihit);
0020
0021
0022 if (reco_lyr < 0)
0023 continue;
0024
0025
0026 if ((reco_idx < 0) || (static_cast<size_t>(reco_idx) >= layerHits[reco_lyr].size()))
0027 continue;
0028
0029
0030 for (int seed_ihit = 0; seed_ihit < seed_trk.nTotalHits(); ++seed_ihit) {
0031 const int seed_lyr = seed_trk.getHitLyr(seed_ihit);
0032 const int seed_idx = seed_trk.getHitIdx(seed_ihit);
0033
0034
0035 if (seed_lyr < 0)
0036 continue;
0037
0038
0039 if (reco_lyr != seed_lyr)
0040 continue;
0041
0042
0043 if ((seed_idx < 0) || (static_cast<size_t>(seed_idx) >= layerHits[seed_lyr].size()))
0044 continue;
0045
0046
0047 if (reco_idx == seed_idx)
0048 matchedSeedHits_.emplace_back(seed_idx, seed_lyr);
0049 }
0050 }
0051 }
0052
0053 bool TrackExtra::isSeedHit(const int lyr, const int idx) const {
0054 return (std::find_if(matchedSeedHits_.begin(), matchedSeedHits_.end(), [=](const auto& matchedSeedHit) {
0055 return ((matchedSeedHit.layer == lyr) && (matchedSeedHit.index == idx));
0056 }) != matchedSeedHits_.end());
0057 }
0058
0059 int TrackExtra::modifyRefTrackID(const int foundHits,
0060 const int minHits,
0061 const TrackVec& reftracks,
0062 const int trueID,
0063 const int duplicate,
0064 int refTrackID) {
0065
0066 if (duplicate) {
0067 refTrackID = -10;
0068 } else {
0069 if (refTrackID >= 0) {
0070 if (reftracks[refTrackID].isFindable()) {
0071 if (foundHits < minHits)
0072 refTrackID = -2;
0073
0074 } else
0075 {
0076 if (foundHits < minHits)
0077 refTrackID = -3;
0078 else
0079 refTrackID = -4;
0080 }
0081 } else if (refTrackID == -1) {
0082 if (trueID >= 0) {
0083 if (reftracks[trueID].isFindable()) {
0084 if (foundHits < minHits)
0085 refTrackID = -5;
0086
0087 } else
0088 {
0089 if (foundHits < minHits)
0090 refTrackID = -6;
0091 else
0092 refTrackID = -7;
0093 }
0094 } else {
0095 if (foundHits < minHits)
0096 refTrackID = -8;
0097 else
0098 refTrackID = -9;
0099 }
0100 }
0101 }
0102 return refTrackID;
0103 }
0104
0105
0106 void TrackExtra::setMCTrackIDInfo(const Track& trk,
0107 const std::vector<HitVec>& layerHits,
0108 const MCHitInfoVec& globalHitInfo,
0109 const TrackVec& simtracks,
0110 const bool isSeed,
0111 const bool isPure) {
0112 dprintf("TrackExtra::setMCTrackIDInfo for track with label %d, total hits %d, found hits %d\n",
0113 trk.label(),
0114 trk.nTotalHits(),
0115 trk.nFoundHits());
0116
0117 std::vector<int> mcTrackIDs;
0118 int nSeedHits = nMatchedSeedHits();
0119
0120
0121 for (int ihit = 0; ihit < trk.nTotalHits(); ++ihit) {
0122 const int lyr = trk.getHitLyr(ihit);
0123 const int idx = trk.getHitIdx(ihit);
0124
0125
0126 if (lyr < 0)
0127 continue;
0128
0129
0130 if (!Config::mtvLikeValidation && !isSeed && isSeedHit(lyr, idx))
0131 continue;
0132
0133
0134 if ((idx >= 0) && (static_cast<size_t>(idx) < layerHits[lyr].size())) {
0135
0136 const int mchitid = layerHits[lyr][idx].mcHitID();
0137 mcTrackIDs.push_back(globalHitInfo[mchitid].mcTrackID());
0138
0139 dprintf(" ihit=%3d trk.hit_idx=%4d trk.hit_lyr=%2d mchitid=%4d mctrkid=%3d\n",
0140 ihit,
0141 idx,
0142 lyr,
0143 mchitid,
0144 globalHitInfo[mchitid].mcTrackID());
0145 } else {
0146 dprintf(" ihit=%3d trk.hit_idx=%4d trk.hit_lyr=%2d\n", ihit, idx, lyr);
0147 }
0148 }
0149
0150 int mccount = 0;
0151 int mcTrackID = -1;
0152 if (!mcTrackIDs.empty())
0153 {
0154
0155 std::sort(mcTrackIDs.begin(), mcTrackIDs.end());
0156
0157
0158 mcTrackIDs.erase(std::remove_if(mcTrackIDs.begin(), mcTrackIDs.end(), [](const int id) { return id < 0; }),
0159 mcTrackIDs.end());
0160
0161 int n_ids = mcTrackIDs.size();
0162 int i = 0;
0163 while (i < n_ids) {
0164 int j = i + 1;
0165 while (j < n_ids && mcTrackIDs[j] == mcTrackIDs[i])
0166 ++j;
0167
0168 int n = j - i;
0169 if (mcTrackIDs[i] >= 0 && n > mccount) {
0170 mcTrackID = mcTrackIDs[i];
0171 mccount = n;
0172 }
0173 i = j;
0174 }
0175
0176
0177 const int nCandHits = ((Config::mtvLikeValidation || isSeed) ? trk.nFoundHits() : trk.nFoundHits() - nSeedHits);
0178
0179
0180 if ((Config::mtvLikeValidation ? (4 * mccount > 3 * nCandHits) : (2 * mccount >= nCandHits))) {
0181
0182 if (isPure) {
0183 if (mcTrackID == seedID_)
0184 mcTrackID_ = mcTrackID;
0185 else
0186 mcTrackID_ = -1;
0187 } else {
0188 mcTrackID_ = mcTrackID;
0189 }
0190 } else
0191 {
0192 mcTrackID_ = -1;
0193 }
0194
0195
0196 if (isPure && mcTrackID == -1) {
0197 mccount = 0;
0198 for (auto id : mcTrackIDs) {
0199 if (id == seedID_)
0200 mccount++;
0201 }
0202 }
0203
0204
0205 nHitsMatched_ = mccount;
0206 fracHitsMatched_ = float(nHitsMatched_) / float(nCandHits);
0207
0208
0209 dPhi_ =
0210 (mcTrackID >= 0 ? squashPhiGeneral(simtracks[mcTrackID].swimPhiToR(trk.x(), trk.y()) - trk.momPhi()) : -99.f);
0211 } else {
0212 mcTrackID_ = mcTrackID;
0213 nHitsMatched_ = -99;
0214 fracHitsMatched_ = -99.f;
0215 dPhi_ = -99.f;
0216 }
0217
0218
0219 if (!isSeed) {
0220 mcTrackID_ = modifyRefTrackID(trk.nFoundHits() - nSeedHits,
0221 Config::nMinFoundHits - nSeedHits,
0222 simtracks,
0223 (isPure ? seedID_ : -1),
0224 trk.getDuplicateValue(),
0225 mcTrackID_);
0226 }
0227
0228 dprint("Track " << trk.label() << " best mc track " << mcTrackID_ << " count " << mccount << "/"
0229 << trk.nFoundHits());
0230 }
0231
0232 typedef std::pair<int, float> idchi2Pair;
0233 typedef std::vector<idchi2Pair> idchi2PairVec;
0234
0235 inline bool sortIDsByChi2(const idchi2Pair& cand1, const idchi2Pair& cand2) { return cand1.second < cand2.second; }
0236
0237 inline int getMatchBin(const float pt) {
0238 if (pt < 0.75f)
0239 return 0;
0240 else if (pt < 1.f)
0241 return 1;
0242 else if (pt < 2.f)
0243 return 2;
0244 else if (pt < 5.f)
0245 return 3;
0246 else if (pt < 10.f)
0247 return 4;
0248 else
0249 return 5;
0250 }
0251
0252 void TrackExtra::setCMSSWTrackIDInfoByTrkParams(const Track& trk,
0253 const std::vector<HitVec>& layerHits,
0254 const TrackVec& cmsswtracks,
0255 const RedTrackVec& redcmsswtracks,
0256 const bool isBkFit) {
0257
0258 const SVector6& trkParams = trk.parameters();
0259 const SMatrixSym66& trkErrs = trk.errors();
0260
0261
0262 const int bin = getMatchBin(trk.pT());
0263
0264
0265 SVector2 trkParamsR;
0266 trkParamsR[0] = trkParams[3];
0267 trkParamsR[1] = trkParams[5];
0268
0269 SMatrixSym22 trkErrsR;
0270 trkErrsR[0][0] = trkErrs[3][3];
0271 trkErrsR[1][1] = trkErrs[5][5];
0272 trkErrsR[0][1] = trkErrs[3][5];
0273 trkErrsR[1][0] = trkErrs[5][3];
0274
0275
0276 idchi2PairVec cands;
0277
0278
0279 for (const auto& redcmsswtrack : redcmsswtracks) {
0280 const float chi2 = std::abs(computeHelixChi2(redcmsswtrack.parameters(), trkParamsR, trkErrsR, false));
0281 if (chi2 < Config::minCMSSWMatchChi2[bin])
0282 cands.push_back(std::make_pair(redcmsswtrack.label(), chi2));
0283 }
0284
0285
0286 float minchi2 = -1e6;
0287 if (!cands.empty()) {
0288 std::sort(cands.begin(), cands.end(), sortIDsByChi2);
0289 minchi2 = cands.front().second;
0290 }
0291
0292
0293 int cmsswTrackID = -1;
0294 int nHitsMatched = 0;
0295 float bestdPhi = Config::minCMSSWMatchdPhi[bin];
0296 float bestchi2 = minchi2;
0297
0298
0299 for (auto&& cand : cands) {
0300
0301 const auto label = cand.first;
0302 const auto& cmsswtrack = cmsswtracks[label];
0303
0304
0305 const float diffPhi =
0306 squashPhiGeneral((isBkFit ? cmsswtrack.momPhi() : cmsswtrack.swimPhiToR(trk.x(), trk.y())) - trk.momPhi());
0307
0308
0309 if (std::abs(diffPhi) < std::abs(bestdPhi)) {
0310 const HitLayerMap& hitLayerMap = redcmsswtracks[label].hitLayerMap();
0311 int matched = 0;
0312
0313
0314 for (int ihit = 0; ihit < trk.nTotalHits(); ihit++) {
0315 const int lyr = trk.getHitLyr(ihit);
0316 const int idx = trk.getHitIdx(ihit);
0317
0318
0319 if (isSeedHit(lyr, idx))
0320 continue;
0321
0322
0323 if (idx < 0 || !hitLayerMap.count(lyr))
0324 continue;
0325
0326
0327 for (auto cidx : hitLayerMap.at(lyr)) {
0328
0329 if (cidx == idx) {
0330 matched++;
0331 break;
0332 }
0333 }
0334 }
0335
0336
0337 bestdPhi = diffPhi;
0338 nHitsMatched = matched;
0339 cmsswTrackID = label;
0340 bestchi2 = cand.second;
0341 }
0342 }
0343
0344
0345 cmsswTrackID_ = cmsswTrackID;
0346 helixChi2_ = bestchi2;
0347 dPhi_ = bestdPhi;
0348
0349
0350 const int nSeedHits = nMatchedSeedHits();
0351
0352
0353 cmsswTrackID_ = modifyRefTrackID(trk.nFoundHits() - nSeedHits,
0354 Config::nMinFoundHits - nSeedHits,
0355 cmsswtracks,
0356 -1,
0357 trk.getDuplicateValue(),
0358 cmsswTrackID_);
0359
0360
0361 nHitsMatched_ = nHitsMatched;
0362 fracHitsMatched_ =
0363 float(nHitsMatched_) / float(trk.nFoundHits() - nSeedHits);
0364 }
0365
0366 void TrackExtra::setCMSSWTrackIDInfoByHits(const Track& trk,
0367 const LayIdxIDVecMapMap& cmsswHitIDMap,
0368 const TrackVec& cmsswtracks,
0369 const TrackExtraVec& cmsswextras,
0370 const RedTrackVec& redcmsswtracks,
0371 const int cmsswlabel) {
0372
0373
0374
0375 std::unordered_map<int, int> labelMatchMap;
0376
0377
0378 for (int ihit = 0; ihit < trk.nTotalHits(); ihit++) {
0379 const int lyr = trk.getHitLyr(ihit);
0380 const int idx = trk.getHitIdx(ihit);
0381
0382 if (lyr < 0 || idx < 0)
0383 continue;
0384 if (isSeedHit(lyr, idx))
0385 continue;
0386 if (!cmsswHitIDMap.count(lyr))
0387 continue;
0388 if (!cmsswHitIDMap.at(lyr).count(idx))
0389 continue;
0390 {
0391 for (const auto label : cmsswHitIDMap.at(lyr).at(idx)) {
0392 labelMatchMap[label]++;
0393 }
0394 }
0395 }
0396
0397
0398 std::vector<int> labelMatchVec;
0399 for (const auto labelMatchPair : labelMatchMap) {
0400 const auto cmsswlabel = labelMatchPair.first;
0401 const auto nMatchedHits = labelMatchPair.second;
0402
0403
0404 if ((2 * nMatchedHits) >= (cmsswtracks[cmsswlabel].nUniqueLayers() - cmsswextras[cmsswlabel].nMatchedSeedHits()))
0405 labelMatchVec.push_back(cmsswlabel);
0406 }
0407
0408
0409 int cmsswTrackID = -1;
0410
0411
0412 if (!labelMatchVec.empty()) {
0413
0414 std::sort(labelMatchVec.begin(), labelMatchVec.end(), [&](const int label1, const int label2) {
0415 if (labelMatchMap[label1] == labelMatchMap[label2]) {
0416 const auto& track1 = cmsswtracks[label1];
0417 const auto& track2 = cmsswtracks[label2];
0418
0419 const auto& extra1 = cmsswextras[label1];
0420 const auto& extra2 = cmsswextras[label2];
0421
0422 return ((track1.nUniqueLayers() - extra1.nMatchedSeedHits()) <
0423 (track2.nUniqueLayers() - extra2.nMatchedSeedHits()));
0424 }
0425 return labelMatchMap[label1] > labelMatchMap[label2];
0426 });
0427
0428
0429 cmsswTrackID = labelMatchVec.front();
0430
0431
0432 if (cmsswlabel >= 0) {
0433 if (cmsswTrackID == cmsswlabel) {
0434 cmsswTrackID_ = cmsswTrackID;
0435 } else {
0436 cmsswTrackID = cmsswlabel;
0437 cmsswTrackID_ = -1;
0438 }
0439 } else
0440 {
0441 cmsswTrackID_ = cmsswTrackID;
0442 }
0443
0444
0445 nHitsMatched_ = labelMatchMap[cmsswTrackID];
0446 } else
0447 {
0448
0449 cmsswTrackID_ = cmsswTrackID;
0450
0451
0452 int nHitsMatched = 0;
0453
0454
0455 if (cmsswlabel >= 0) {
0456 cmsswTrackID = cmsswlabel;
0457 nHitsMatched = labelMatchMap[cmsswTrackID];
0458 } else {
0459
0460 for (const auto labelMatchPair : labelMatchMap) {
0461 if (labelMatchPair.second > nHitsMatched) {
0462 cmsswTrackID = labelMatchPair.first;
0463 nHitsMatched = labelMatchPair.second;
0464 }
0465 }
0466 }
0467
0468 nHitsMatched_ = nHitsMatched;
0469 }
0470
0471
0472 if (cmsswTrackID >= 0) {
0473
0474 const SVector6& trkParams = trk.parameters();
0475 const SMatrixSym66& trkErrs = trk.errors();
0476
0477
0478 SVector2 trkParamsR;
0479 trkParamsR[0] = trkParams[3];
0480 trkParamsR[1] = trkParams[5];
0481
0482 SMatrixSym22 trkErrsR;
0483 trkErrsR[0][0] = trkErrs[3][3];
0484 trkErrsR[1][1] = trkErrs[5][5];
0485 trkErrsR[0][1] = trkErrs[3][5];
0486 trkErrsR[1][0] = trkErrs[5][3];
0487
0488
0489 helixChi2_ = std::abs(computeHelixChi2(redcmsswtracks[cmsswTrackID].parameters(), trkParamsR, trkErrsR, false));
0490 dPhi_ = squashPhiGeneral(cmsswtracks[cmsswTrackID].swimPhiToR(trk.x(), trk.y()) - trk.momPhi());
0491 } else {
0492 helixChi2_ = -99.f;
0493 dPhi_ = -99.f;
0494 }
0495
0496
0497 const int nSeedHits = nMatchedSeedHits();
0498
0499
0500 cmsswTrackID_ = modifyRefTrackID(trk.nFoundHits() - nSeedHits,
0501 Config::nMinFoundHits - nSeedHits,
0502 cmsswtracks,
0503 cmsswlabel,
0504 trk.getDuplicateValue(),
0505 cmsswTrackID_);
0506
0507
0508 fracHitsMatched_ = (cmsswTrackID >= 0 ? (float(nHitsMatched_) / float(cmsswtracks[cmsswTrackID].nUniqueLayers() -
0509 cmsswextras[cmsswTrackID].nMatchedSeedHits()))
0510 : 0.f);
0511 }
0512
0513 }