PixelTrackFitting/src/PixelTrackCleanerBySharedHits.cc

0001 #include "RecoTracker/PixelTrackFitting/interface/PixelTrackCleanerBySharedHits.h"
0002
0003 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0004 #include "DataFormats/TrackReco/interface/Track.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0007 #include <cassert>
0008 #include <cstdint>
0009
0010 using namespace std;
0011 using namespace reco;
0012 using namespace pixeltrackfitting;
0013
0014 PixelTrackCleanerBySharedHits::PixelTrackCleanerBySharedHits(bool useQuadrupletAlgo)
0015     : PixelTrackCleaner(true),  // to mark this as fast algo
0016       useQuadrupletAlgo_(useQuadrupletAlgo) {}
0017
0018 PixelTrackCleanerBySharedHits::~PixelTrackCleanerBySharedHits() {}
0019
0020 void PixelTrackCleanerBySharedHits::cleanTracks(TracksWithTTRHs& trackHitPairs) const {
0021   LogDebug("PixelTrackCleanerBySharedHits") << "Cleanering tracks"
0022                                             << "\n";
0023   unsigned int size = trackHitPairs.size();
0024   if (size <= 1)
0025     return;
0026
0027   // sort (stabilize cleaning)
0028   uint16_t score[size];
0029   unsigned int ind[size];
0030   for (auto i = 0U; i < size; ++i) {
0031     ind[i] = i;
0032     score[i] = 32000 - std::min(32000, int(trackHitPairs[i].first->chi2() * 100.f));  // chi2: smaller is better
0033     if (trackHitPairs[i].second.size() == 4)
0034       score[i] += 32001;  //  s4 always better than s3
0035   }
0036   std::sort(ind, ind + size, [&](unsigned int i, unsigned int j) { return score[i] > score[j]; });
0037
0038   auto kill = [&](unsigned int k) {
0039     delete trackHitPairs[k].first;
0040     trackHitPairs[k].first = nullptr;
0041   };
0042
0043   // sorted: first is always better!
0044
0045   // first loop: only first    two hits....
0046   for (auto i = 0U; i < size; ++i) {
0047     auto iTrack1 = ind[i];
0048     auto track1 = trackHitPairs[iTrack1].first;
0049     if (!track1)
0050       continue;
0051     auto const& recHits1 = trackHitPairs[iTrack1].second;
0052     for (auto j = i + 1; j < size; ++j) {
0053       auto iTrack2 = ind[j];
0054       auto track2 = trackHitPairs[iTrack2].first;
0055       if (!track2)
0056         continue;
0057       auto const& recHits2 = trackHitPairs[iTrack2].second;
0058       if (recHits1[0] != recHits2[0])
0059         continue;
0060       if (recHits1[1] != recHits2[1])
0061         continue;
0062       kill(iTrack2);
0063     }  // tk2
0064   }    // tk1
0065
0066   // second loop: first and third hits....
0067   for (auto i = 0U; i < size; ++i) {
0068     auto iTrack1 = ind[i];
0069     auto track1 = trackHitPairs[iTrack1].first;
0070     if (!track1)
0071       continue;
0072     auto const& recHits1 = trackHitPairs[iTrack1].second;
0073     if (recHits1.size() < 3)
0074       continue;
0075     for (auto j = i + 1; j < size; ++j) {
0076       auto iTrack2 = ind[j];
0077       auto track2 = trackHitPairs[iTrack2].first;
0078       if (!track2)
0079         continue;
0080       auto const& recHits2 = trackHitPairs[iTrack2].second;
0081       if (recHits2.size() < 3)
0082         continue;
0083       if (recHits1[0] != recHits2[0])
0084         continue;
0085       if (recHits1[2] != recHits2[2])
0086         continue;
0087       kill(iTrack2);
0088     }  // tk2
0089   }    // tk1
0090
0091   // final loop: all the rest
0092   for (auto i = 0U; i < size; ++i) {
0093     auto iTrack1 = ind[i];
0094     auto track1 = trackHitPairs[iTrack1].first;
0095     if (!track1)
0096       continue;
0097     auto const& recHits1 = trackHitPairs[iTrack1].second;
0098     auto s1 = recHits1.size();
0099     for (auto j = i + 1; j < size; ++j) {
0100       auto iTrack2 = ind[j];
0101       auto track2 = trackHitPairs[iTrack2].first;
0102       if (!track2)
0103         continue;
0104       auto const& recHits2 = trackHitPairs[iTrack2].second;
0105       auto s2 = recHits2.size();
0106       auto commonRecHits = 0U;
0107       auto f2 = 0U;
0108       for (auto iRecHit1 = 0U; iRecHit1 < s1; ++iRecHit1) {
0109         for (auto iRecHit2 = f2; iRecHit2 < s2; ++iRecHit2) {
0110           if (recHits1[iRecHit1] == recHits2[iRecHit2]) {
0111             ++commonRecHits;
0112             f2 = iRecHit2 + 1;
0113             break;
0114           }  // if a hit is common, no other can be the same!
0115         }
0116         if (commonRecHits > 1)
0117           break;
0118       }
0119       if (useQuadrupletAlgo_) {
0120         if (commonRecHits >= 1) {
0121           if (s1 == 3 || commonRecHits > 1) {
0122             kill(iTrack2);
0123           }
0124         }
0125       } else if (commonRecHits > 1) {
0126         kill(iTrack2);
0127       }
0128     }  // tk2
0129   }    //tk1
0130
0131   trackHitPairs.erase(
0132       std::remove_if(trackHitPairs.begin(), trackHitPairs.end(), [&](TrackWithTTRHs& v) { return nullptr == v.first; }),
0133       trackHitPairs.end());
0134 }