PixelLowPtUtilities/src/TrackCleaner.cc

0001 #include "RecoTracker/PixelLowPtUtilities/interface/TrackCleaner.h"
0002 #include "RecoTracker/PixelLowPtUtilities/interface/HitInfo.h"
0003
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0007 #include "DataFormats/TrackReco/interface/Track.h"
0008
0009 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0010 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0011 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0012
0013 using namespace std;
0014 using namespace pixeltrackfitting;
0015
0016 /*****************************************************************************/
0017 class HitComparatorByRadius {  // No access to geometry!
0018 public:
0019   HitComparatorByRadius(const TrackerTopology *tTopo) { tTopo_ = tTopo; }
0020
0021 private:
0022   const TrackerTopology *tTopo_;
0023
0024 public:
0025   bool operator()(const TrackingRecHit *a, const TrackingRecHit *b) const {
0026     DetId i1 = a->geographicalId();
0027     DetId i2 = b->geographicalId();
0028
0029     bool isBarrel = (i1.subdetId() == int(PixelSubdetector::PixelBarrel));
0030
0031     if (i1.subdetId() != i2.subdetId()) {
0032       return isBarrel;
0033     } else {
0034       if (isBarrel) {
0035         int r1 = (tTopo_->pxbLayer(i1) - 1) * 2 + (tTopo_->pxbLadder(i1) - 1) % 2;
0036         int r2 = (tTopo_->pxbLayer(i2) - 1) * 2 + (tTopo_->pxbLadder(i2) - 1) % 2;
0037
0038         return (r1 < r2);
0039       } else {
0040         int r1 = (tTopo_->pxfDisk(i1) - 1) * 2 + (tTopo_->pxfPanel(i1) - 1) % 2;
0041         int r2 = (tTopo_->pxfDisk(i2) - 1) * 2 + (tTopo_->pxfPanel(i2) - 1) % 2;
0042
0043         return (r1 < r2);
0044       }
0045     }
0046   }
0047 };
0048
0049 /*****************************************************************************/
0050 class HitComparator {
0051 public:
0052   bool operator()(const TrackingRecHit *a, const TrackingRecHit *b) const {
0053     if (a->geographicalId() < b->geographicalId())
0054       return true;
0055     if (b->geographicalId() < a->geographicalId())
0056       return false;
0057
0058     if (a->localPosition().x() < b->localPosition().x() - 1e-5)
0059       return true;
0060     if (b->localPosition().x() < a->localPosition().x() - 1e-5)
0061       return false;
0062
0063     if (a->localPosition().y() < b->localPosition().y() - 1e-5)
0064       return true;
0065     if (b->localPosition().y() < a->localPosition().y() - 1e-5)
0066       return false;
0067
0068     return false;
0069   }
0070 };
0071
0072 /*****************************************************************************/
0073 TrackCleaner::TrackCleaner(const TrackerTopology *tTopo) : tTopo_(tTopo) {}
0074
0075 /*****************************************************************************/
0076 TrackCleaner::~TrackCleaner() {}
0077
0078 /*****************************************************************************/
0079 bool TrackCleaner::areSame(const TrackingRecHit *a, const TrackingRecHit *b) const {
0080   if (a->geographicalId() != b->geographicalId())
0081     return false;
0082
0083   if (fabs(a->localPosition().x() - b->localPosition().x()) < 1e-5 &&
0084       fabs(a->localPosition().y() - b->localPosition().y()) < 1e-5)
0085     return true;
0086   else
0087     return false;
0088 }
0089
0090 /*****************************************************************************/
0091 bool TrackCleaner::isCompatible(const DetId &i1, const DetId &i2) const {
0092   // different subdet
0093   if (i1.subdetId() != i2.subdetId())
0094     return true;
0095
0096   if (i1.subdetId() == int(PixelSubdetector::PixelBarrel)) {  // barrel
0097
0098     if (tTopo_->pxbLayer(i1) != tTopo_->pxbLayer(i2))
0099       return true;
0100
0101     int dphi = abs(int(tTopo_->pxbLadder(i1) - tTopo_->pxbLadder(i2)));
0102     //FIXME: non-phase-0 wrap-around is wrong (sh/be 12, 28, 44, 64 in phase-1)
0103     // the harm is somewhat minimal with some excess of duplicates.
0104     constexpr int max[4] = {20, 32, 44, 64};
0105     auto aLayer = tTopo_->pxbLayer(i1) - 1;
0106     assert(aLayer < 4);
0107     if (dphi > max[aLayer] / 2)
0108       dphi = max[aLayer] - dphi;
0109
0110     int dz = abs(int(tTopo_->pxbModule(i1) - tTopo_->pxbModule(i2)));
0111
0112     if (dphi == 1 && dz <= 1)
0113       return true;
0114   } else {  // endcap
0115
0116     if (tTopo_->pxfSide(i1) != tTopo_->pxfSide(i2) || tTopo_->pxfDisk(i1) != tTopo_->pxfDisk(i2))
0117       return true;
0118
0119     int dphi = abs(int(tTopo_->pxfBlade(i1) - tTopo_->pxfBlade(i2)));
0120     constexpr int max = 24;  //FIXME: non-phase-0 wrap-around is wrong
0121     if (dphi > max / 2)
0122       dphi = max - dphi;
0123
0124     int dr = abs(int(((tTopo_->pxfModule(i1) - 1) * 2 + (tTopo_->pxfPanel(i1) - 1)) -
0125                      ((tTopo_->pxfModule(i2) - 1) * 2 + (tTopo_->pxfPanel(i2) - 1))));
0126
0127     if (dphi <= 1 && dr <= 1 && !(dphi == 0 && dr == 0))
0128       return true;
0129   }
0130
0131   return false;
0132 }
0133
0134 /*****************************************************************************/
0135 bool TrackCleaner::canBeMerged(const vector<const TrackingRecHit *> &recHitsA,
0136                                const vector<const TrackingRecHit *> &recHitsB) const {
0137   bool ok = true;
0138
0139   for (vector<const TrackingRecHit *>::const_iterator recHitA = recHitsA.begin(); recHitA != recHitsA.end(); recHitA++)
0140     for (vector<const TrackingRecHit *>::const_iterator recHitB = recHitsB.begin(); recHitB != recHitsB.end();
0141          recHitB++)
0142       if (!areSame(*recHitA, *recHitB))
0143         if (!isCompatible((*recHitA)->geographicalId(), (*recHitB)->geographicalId()))
0144           ok = false;
0145
0146   return ok;
0147 }
0148
0149 /*****************************************************************************/
0150 TracksWithRecHits TrackCleaner::cleanTracks(const TracksWithRecHits &tracks_) const {
0151   // Local copy
0152   TracksWithRecHits tracks = tracks_;
0153
0154   typedef map<const TrackingRecHit *, vector<unsigned int>, HitComparator> RecHitMap;
0155
0156   vector<bool> keep(tracks.size(), true);
0157
0158   int changes;
0159
0160   LogTrace("MinBiasTracking") << " [TrackCleaner] initial tracks : " << tracks.size();
0161
0162   for (unsigned int i = 0; i < tracks.size(); i++)
0163     LogTrace("TrackCleaner") << "   Track #" << i << " : " << HitInfo::getInfo(tracks[i].second, tTopo_);
0164
0165   do {
0166     changes = 0;
0167
0168     RecHitMap recHitMap;
0169
0170     /*
0171   LogTrace("MinBiasTracking")
0172     << " [TrackCleaner] fill rechit map";
0173 */
0174
0175     // Fill the rechit map
0176     for (unsigned int i = 0; i < tracks.size(); i++)
0177       if (keep[i]) {
0178         for (vector<const TrackingRecHit *>::const_iterator recHit = tracks[i].second.begin();
0179              recHit != tracks[i].second.end();
0180              recHit++)
0181           recHitMap[*recHit].push_back(i);
0182       }
0183
0184     // Look at each track
0185     typedef map<unsigned int, int, less<unsigned int> > TrackMap;
0186
0187     for (unsigned int i = 0; i < tracks.size(); i++) {
0188       // Skip if 'i' already removed
0189       if (!keep[i])
0190         continue;
0191
0192       /*
0193
0194     bool addedNewHit = false;
0195     do
0196     {
0197 */
0198       TrackMap trackMap;
0199
0200       // Go trough all rechits of this track
0201       for (vector<const TrackingRecHit *>::const_iterator recHit = tracks[i].second.begin();
0202            recHit != tracks[i].second.end();
0203            recHit++) {
0204         // Get tracks sharing this rechit
0205         vector<unsigned int> sharing = recHitMap[*recHit];
0206
0207         for (vector<unsigned int>::iterator j = sharing.begin(); j != sharing.end(); j++)
0208           if (i < *j)
0209             trackMap[*j]++;
0210       }
0211
0212       // Check for tracks with shared rechits
0213       for (TrackMap::iterator sharing = trackMap.begin(); sharing != trackMap.end(); sharing++) {
0214         unsigned int j = (*sharing).first;
0215         if (!keep[i] || !keep[j])
0216           continue;
0217
0218         if (tracks[i].second.size() >= 3) {  // triplet tracks
0219           if ((*sharing).second > min(int(tracks[i].second.size()),
0220                                       int(tracks[j].second.size())) /
0221                                       2) {                          // more than min(hits1,hits2)/2 rechits are shared
0222             if (canBeMerged(tracks[i].second, tracks[j].second)) {  // no common layer
0223               // merge tracks, add separate hits of the second to the first one
0224               for (vector<const TrackingRecHit *>::const_iterator recHit = tracks[j].second.begin();
0225                    recHit != tracks[j].second.end();
0226                    recHit++) {
0227                 bool ok = true;
0228                 for (vector<const TrackingRecHit *>::const_iterator recHitA = tracks[i].second.begin();
0229                      recHitA != tracks[i].second.end();
0230                      recHitA++)
0231                   if (areSame(*recHit, *recHitA))
0232                     ok = false;
0233
0234                 if (ok) {
0235                   tracks[i].second.push_back(*recHit);
0236                   recHitMap[*recHit].push_back(i);
0237
0238                   sort(tracks[i].second.begin(), tracks[i].second.end(), HitComparatorByRadius(tTopo_));
0239
0240                   //addedNewHit = true;
0241                 }
0242               }
0243
0244               LogTrace("TrackCleaner") << "   Merge #" << i << " #" << j << ", first now has "
0245                                        << tracks[i].second.size();
0246
0247               // Remove second track
0248               keep[j] = false;
0249
0250               changes++;
0251             } else {  // there is a common layer, keep smaller impact
0252               if (fabs(tracks[i].first->d0()) < fabs(tracks[j].first->d0()))
0253                 keep[j] = false;
0254               else
0255                 keep[i] = false;
0256
0257               LogTrace("TrackCleaner") << "   Clash #" << i << " #" << j << " keep lower d0 " << tracks[i].first->d0()
0258                                        << " " << tracks[j].first->d0() << ", keep #"
0259                                        << (keep[i] ? i : (keep[j] ? j : 9999));
0260
0261               changes++;
0262             }
0263           } else {  // note more than 50%, but at least two are shared
0264             if ((*sharing).second > 1) {
0265               if (tracks[i].second.size() != tracks[j].second.size()) {  // keep longer
0266                 if (tracks[i].second.size() > tracks[j].second.size())
0267                   keep[j] = false;
0268                 else
0269                   keep[i] = false;
0270                 changes++;
0271
0272                 LogTrace("TrackCleaner") << "   Sharing " << (*sharing).second << " remove by size";
0273               } else {  // keep smaller impact
0274                 if (fabs(tracks[i].first->d0()) < fabs(tracks[j].first->d0()))
0275                   keep[j] = false;
0276                 else
0277                   keep[i] = false;
0278                 changes++;
0279
0280                 LogTrace("TrackCleaner") << "   Sharing " << (*sharing).second << " remove by d0";
0281               }
0282             }
0283           }
0284         } else {  // pair tracks
0285           if ((*sharing).second > 0) {
0286             // Remove second track
0287             keep[j] = false;
0288
0289             changes++;
0290           }
0291         }
0292       }
0293       /*
0294     }
0295     while(addedNewHit);
0296 */
0297     }
0298   } while (changes > 0);
0299
0300   // Final copy
0301   TracksWithRecHits cleaned;
0302
0303   for (unsigned int i = 0; i < tracks.size(); i++)
0304     if (keep[i])
0305       cleaned.push_back(tracks[i]);
0306     else
0307       delete tracks_[i].first;
0308
0309   LogTrace("MinBiasTracking") << " [TrackCleaner] cleaned tracks : " << cleaned.size();
0310
0311   for (unsigned int i = 0; i < cleaned.size(); i++)
0312     LogTrace("TrackCleaner") << "   Track #" << i << " : " << HitInfo::getInfo(cleaned[i].second, tTopo_);
0313
0314   return cleaned;
0315 }