TrackFindingTMTT/src/DupFitTrkKiller.cc

0001 #include "L1Trigger/TrackFindingTMTT/interface/DupFitTrkKiller.h"
0002 #include "L1Trigger/TrackFindingTMTT/interface/Settings.h"
0003 #include "L1Trigger/TrackFindingTMTT/interface/PrintL1trk.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005 #include <map>
0006
0007 using namespace std;
0008
0009 namespace tmtt {
0010
0011   //=== Make available cfg parameters & specify which algorithm is to be used for duplicate track removal.
0012
0013   DupFitTrkKiller::DupFitTrkKiller(const Settings* settings)
0014       : settings_(settings), dupTrkAlg_(static_cast<DupAlgoName>(settings->dupTrkAlgFit())) {}
0015
0016   //=== Eliminate duplicate tracks from the input collection, and so return a reduced list of tracks.
0017
0018   list<const L1fittedTrack*> DupFitTrkKiller::filter(const list<L1fittedTrack>& vecTracks) const {
0019     if (dupTrkAlg_ == DupAlgoName::None) {
0020       // We are not running duplicate removal, so return original fitted track collection.
0021       list<const L1fittedTrack*> copyTracks;
0022       for (const L1fittedTrack& trk : vecTracks) {
0023         copyTracks.push_back(&trk);
0024       }
0025       return copyTracks;
0026
0027     } else {
0028       // Choose which algorithm to run, based on parameter dupTrkAlg_.
0029       switch (dupTrkAlg_) {
0030         // Run filters that only work on fitted tracks.
0031         case DupAlgoName::Algo1:
0032           return filterAlg1(vecTracks);
0033           break;
0034         case DupAlgoName::Algo2:
0035           return filterAlg2(vecTracks);
0036           break;
0037         default:
0038           throw cms::Exception("BadConfig") << "KillDupTrks: Cfg option DupFitTrkAlg has invalid value.";
0039       }
0040     }
0041   }
0042
0043   //=== Duplicate removal algorithm designed to run after the track helix fit, which eliminates duplicates
0044   //=== simply by requiring that the fitted (q/Pt, phi0) of the track correspond to the same HT cell in
0045   //=== which the track was originally found by the HT.
0046   //=== N.B. This code runs on tracks in a single sector. It could be extended to run on tracks in entire
0047   //=== tracker by adding the track's sector number to memory "htCellUsed" below.
0048
0049   list<const L1fittedTrack*> DupFitTrkKiller::filterAlg1(const list<L1fittedTrack>& tracks) const {
0050     // Hard-wired options to play with.
0051     constexpr bool debug = false;
0052     constexpr bool doRecoveryStep = true;  // Do 2nd pass through rejected tracks to see if any should be rescued.
0053     constexpr bool reduceDups = true;      // Option attempting to reduce duplicate tracks during 2nd pass.
0054     constexpr bool memorizeAllHTcells =
0055         false;  // First pass stores in memory all cells that the HT found tracks in, not just those of tracks accepted by the first pass.
0056     constexpr bool doSectorCheck = false;  // Require fitted helix to lie within sector.
0057     constexpr bool usePtAndZ0Cuts = false;
0058     constexpr bool goOutsideArray = true;  // Also store in memory stubs outside the HT array during 2nd pass.
0059     constexpr bool limitDiff = true;       // Limit allowed diff. between HT & Fit cell to <= 1.
0060
0061     if (debug && not tracks.empty())
0062       PrintL1trk() << "Start DupFitTrkKiller" << tracks.size();
0063
0064     list<const L1fittedTrack*> tracksFiltered;
0065
0066     // Make a first pass through the tracks, doing initial identification of duplicate tracks.
0067     // N.B. BY FILLING THIS WITH CELLS AROUND SELECTED TRACKS, RATHER THAN JUST THE CELL CONTAINING THE
0068     // TRACK, ONE CAN REDUCE THE DUPLICATE RATE FURTHER, AT COST TO EFFICIENCY.
0069     set<pair<unsigned int, unsigned int>> htCellUsed;
0070     list<const L1fittedTrack*> tracksRejected;
0071
0072     // For checking if multiple tracks corresponding to same TP are accepted by duplicate removal.
0073     map<unsigned int, pair<unsigned int, unsigned int>> tpFound;
0074     map<unsigned int, unsigned int> tpFoundAtPass;
0075
0076     for (const L1fittedTrack& trk : tracks) {
0077       // Only consider tracks whose fitted helix parameters are in the same sector as the HT originally used to find the track.
0078       if ((!doSectorCheck) || trk.consistentSector()) {
0079         if ((!usePtAndZ0Cuts) ||
0080             (std::abs(trk.z0()) < settings_->beamWindowZ() && trk.pt() > settings_->houghMinPt() - 0.2)) {
0081           // For debugging.
0082           const TP* tp = trk.matchedTP();
0083
0084           // Check if this track's fitted (q/pt, phi0) helix parameters correspond to the same HT cell as the HT originally found the track in.
0085           bool consistentCell = trk.consistentHTcell();
0086           if (consistentCell) {
0087             // This track is probably not a duplicate, so keep & and store its HT cell location (which equals the HT cell corresponding to the fitted track).
0088             tracksFiltered.push_back(&trk);
0089             // Memorize HT cell location corresponding to this track (identical for HT track & fitted track).
0090             if (!memorizeAllHTcells) {
0091               pair<unsigned int, unsigned int> htCell = trk.cellLocationHT();
0092               htCellUsed.insert(htCell);
0093               if (trk.l1track3D()->mergedHTcell()) {
0094                 // If this is a merged cell, block the other elements too, in case a track found by the HT in an unmerged cell
0095                 // has a fitted cell there.
0096                 pair<unsigned int, unsigned int> htCell10(htCell.first + 1, htCell.second);
0097                 pair<unsigned int, unsigned int> htCell01(htCell.first, htCell.second + 1);
0098                 pair<unsigned int, unsigned int> htCell11(htCell.first + 1, htCell.second + 1);
0099                 htCellUsed.insert(htCell10);
0100                 htCellUsed.insert(htCell01);
0101                 htCellUsed.insert(htCell11);
0102               }
0103             }
0104
0105             if (debug && tp != nullptr) {
0106               PrintL1trk() << "FIRST PASS: m=" << trk.cellLocationHT().first << "/" << trk.cellLocationFit().first
0107                            << " c=" << trk.cellLocationHT().second << "/" << trk.cellLocationFit().second
0108                            << " Delta(m,c)=(" << int(trk.cellLocationHT().first) - int(trk.cellLocationFit().first)
0109                            << "," << int(trk.cellLocationHT().second) - int(trk.cellLocationFit().second)
0110                            << ") pure=" << trk.purity() << " merged=" << trk.l1track3D()->mergedHTcell()
0111                            << " #layers=" << trk.l1track3D()->numLayers() << " tp=" << tp->index() << " dupCell=("
0112                            << tpFound[tp->index()].first << "," << tpFound[tp->index()].second
0113                            << ") dup=" << tpFoundAtPass[tp->index()];
0114               // If the following two variables are non-zero in printout, then track has already been found,
0115               // so we have mistakenly kept a duplicate.
0116               if (tpFound.find(tp->index()) != tpFound.end())
0117                 tpFound[tp->index()] = trk.cellLocationFit();
0118               tpFoundAtPass[tp->index()] = 1;
0119             }
0120
0121           } else {
0122             if (limitDiff) {
0123               const unsigned int maxDiff = 1;
0124               if (abs(int(trk.cellLocationHT().first) - int(trk.cellLocationFit().first)) <= int(maxDiff) &&
0125                   abs(int(trk.cellLocationHT().second) - int(trk.cellLocationFit().second)) <= int(maxDiff))
0126                 tracksRejected.push_back(&trk);
0127             } else {
0128               tracksRejected.push_back(&trk);
0129             }
0130
0131             if (debug && tp != nullptr) {
0132               PrintL1trk() << "FIRST REJECT: m=" << trk.cellLocationHT().first << "/" << trk.cellLocationFit().first
0133                            << " c=" << trk.cellLocationHT().second << "/" << trk.cellLocationFit().second
0134                            << " Delta(m,c)=(" << int(trk.cellLocationHT().first) - int(trk.cellLocationFit().first)
0135                            << "," << int(trk.cellLocationHT().second) - int(trk.cellLocationFit().second)
0136                            << ") pure=" << trk.purity() << " merged=" << trk.l1track3D()->mergedHTcell()
0137                            << " #layers=" << trk.l1track3D()->numLayers() << " tp=" << tp->index() << " dupCell=("
0138                            << tpFound[tp->index()].first << "," << tpFound[tp->index()].second
0139                            << ") dup=" << tpFoundAtPass[tp->index()];
0140             }
0141           }
0142           // Memorize HT cell location corresponding to this track, even if it was not accepted by first pass..
0143           if (memorizeAllHTcells) {
0144             pair<unsigned int, unsigned int> htCell =
0145                 trk.cellLocationFit();  // Intentionally used fit instead of HT here.
0146             htCellUsed.insert(htCell);
0147             if (trk.l1track3D()->mergedHTcell()) {
0148               // If this is a merged cell, block the other elements too, in case a track found by the HT in an unmerged cell
0149               // has a fitted cell there.
0150               // N.B. NO GOOD REASON WHY "-1" IS NOT DONE HERE TOO. MIGHT REDUCE DUPLICATE RATE?
0151               pair<unsigned int, unsigned int> htCell10(htCell.first + 1, htCell.second);
0152               pair<unsigned int, unsigned int> htCell01(htCell.first, htCell.second + 1);
0153               pair<unsigned int, unsigned int> htCell11(htCell.first + 1, htCell.second + 1);
0154               htCellUsed.insert(htCell10);
0155               htCellUsed.insert(htCell01);
0156               htCellUsed.insert(htCell11);
0157             }
0158           }
0159         }
0160       }
0161     }
0162
0163     if (doRecoveryStep) {
0164       // Making a second pass through the rejected tracks, checking if any should be rescued.
0165       for (const L1fittedTrack* trk : tracksRejected) {
0166         // Get location in HT array corresponding to fitted track helix parameters.
0167         pair<unsigned int, unsigned int> htCell = trk->cellLocationFit();
0168         // If this HT cell was not already memorized, rescue this track, since it is probably not a duplicate,
0169         // but just a track whose fitted helix parameters are a bit wierd for some reason.
0170         if (htCellUsed.count(htCell) == 0) {
0171           tracksFiltered.push_back(trk);  // Rescue track.
0172           // Optionally store cell location to avoid rescuing other tracks at the same location, which may be duplicates of this track.
0173           bool outsideCheck = (goOutsideArray || trk->pt() > settings_->houghMinPt());
0174           if (reduceDups && outsideCheck)
0175             htCellUsed.insert(htCell);
0176
0177           // For debugging.
0178           const TP* tp = trk->matchedTP();
0179
0180           if (debug && tp != nullptr) {
0181             PrintL1trk() << "SECOND PASS: m=" << trk->cellLocationHT().first << "/" << trk->cellLocationFit().first
0182                          << " c=" << trk->cellLocationHT().second << "/" << trk->cellLocationFit().second
0183                          << " Delta(m,c)=(" << int(trk->cellLocationHT().first) - int(trk->cellLocationFit().first)
0184                          << "," << int(trk->cellLocationHT().second) - int(trk->cellLocationFit().second)
0185                          << ") pure=" << trk->purity() << " merged=" << trk->l1track3D()->mergedHTcell()
0186                          << " #layers=" << trk->l1track3D()->numLayers() << " tp=" << tp->index() << " dupCell=("
0187                          << tpFound[tp->index()].first << "," << tpFound[tp->index()].second
0188                          << ") dup=" << tpFoundAtPass[tp->index()];
0189             if (tpFound.find(tp->index()) != tpFound.end())
0190               tpFound[tp->index()] = htCell;
0191             tpFoundAtPass[tp->index()] = 2;
0192           }
0193         }
0194       }
0195     }
0196
0197     // Debug printout to identify duplicate tracks that survived.
0198     if (debug)
0199       this->printDuplicateTracks(tracksFiltered);
0200
0201     return tracksFiltered;
0202   }
0203
0204   //=== Duplicate removal algorithm designed to run after the track helix fit, which eliminates duplicates
0205   //=== simply by requiring that no two tracks should have fitted (q/Pt, phi0) that correspond to the same HT
0206   //=== cell. If they do, then only the first to arrive is kept.
0207   //=== N.B. This code runs on tracks in a single sector. It could be extended to run on tracks in entire
0208   //=== tracker by adding the track's sector number to memory "htCellUsed" below.
0209
0210   list<const L1fittedTrack*> DupFitTrkKiller::filterAlg2(const list<L1fittedTrack>& tracks) const {
0211     // Hard-wired options to play with.
0212     const bool debug = false;
0213
0214     if (debug && not tracks.empty())
0215       PrintL1trk() << "START " << tracks.size();
0216
0217     list<const L1fittedTrack*> tracksFiltered;
0218     set<pair<unsigned int, unsigned int>> htCellUsed;
0219
0220     for (const L1fittedTrack& trk : tracks) {
0221       // Get location in HT array corresponding to fitted track helix parameters.
0222       pair<unsigned int, unsigned int> htCell = trk.cellLocationFit();
0223       // If this HT cell was not already memorized, rescue this track, since it is probably not a duplicate,
0224       // but just a track whose fitted helix parameters are a bit wierd for some reason.
0225       if (htCellUsed.count(htCell) == 0) {
0226         tracksFiltered.push_back(&trk);  // Rescue track.
0227         // Store cell location to avoid rescuing other tracks at the same location, which may be duplicates of this track.
0228         htCellUsed.insert(htCell);
0229         if (debug) {
0230           const TP* tp = trk.matchedTP();
0231           int tpIndex = (tp != nullptr) ? tp->index() : -999;
0232           PrintL1trk() << "ALG51: m=" << trk.cellLocationHT().first << "/" << trk.cellLocationFit().first
0233                        << " c=" << trk.cellLocationHT().second << "/" << trk.cellLocationFit().second
0234                        << " tp=" << tpIndex << " pure=" << trk.purity();
0235         }
0236       }
0237     }
0238
0239     // Debug printout to identify duplicate tracks that survived.
0240     if (debug)
0241       this->printDuplicateTracks(tracksFiltered);
0242
0243     return tracksFiltered;
0244   }
0245
0246   // Debug printout of which tracks are duplicates.
0247   void DupFitTrkKiller::printDuplicateTracks(const list<const L1fittedTrack*>& tracks) const {
0248     map<const TP*, list<const L1fittedTrack*>> tpMap;
0249     for (const L1fittedTrack* trk : tracks) {
0250       const TP* tp = trk->matchedTP();
0251       if (tp != nullptr) {
0252         tpMap[tp].push_back(trk);
0253       }
0254     }
0255     for (const auto& p : tpMap) {
0256       const TP* tp = p.first;
0257       const list<const L1fittedTrack*> vecTrk = p.second;
0258       if (vecTrk.size() > 1) {
0259         for (const L1fittedTrack* trk : vecTrk) {
0260           PrintL1trk() << "  MESS UP : m=" << trk->cellLocationHT().first << "/" << trk->cellLocationFit().first
0261                        << " c=" << trk->cellLocationHT().second << "/" << trk->cellLocationFit().second
0262                        << " tp=" << tp->index() << " tp_pt=" << tp->pt() << " fit_pt=" << trk->pt()
0263                        << " pure=" << trk->purity();
0264           PrintL1trk() << "     stubs = ";
0265           for (const Stub* s : trk->stubs())
0266             PrintL1trk() << s->index() << " ";
0267           PrintL1trk();
0268         }
0269       }
0270     }
0271     if (not tracks.empty())
0272       PrintL1trk() << "FOUND " << tracks.size();
0273   }
0274
0275 }  // namespace tmtt