TrajectoryCleaning/src/TrajectoryCleanerBySharedHits.cc

0001 #include "TrackingTools/TrajectoryCleaning/interface/TrajectoryCleanerBySharedHits.h"
0002 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0003 #include "TrackingTools/TransientTrackingRecHit/interface/RecHitComparatorByPosition.h"
0004
0005 #include "TrackingTools/TrajectoryCleaning/src/OtherHashMaps.h"
0006
0007 //#define DEBUG_PRINT(X) X
0008 #define DEBUG_PRINT(X)
0009
0010 namespace {
0011
0012   // Define when two rechits are equals
0013   struct EqualsBySharesInput {
0014     bool operator()(const TransientTrackingRecHit *h1, const TransientTrackingRecHit *h2) const {
0015       return (h1 == h2) || ((h1->geographicalId() == h2->geographicalId()) &&
0016                             (h1->hit()->sharesInput(h2->hit(), TrackingRecHit::some)));
0017     }
0018   };
0019   // Define a hash, i.e. a number that must be equal if hits are equal, and should be different if they're not
0020   struct HashByDetId {
0021     std::size_t operator()(const TransientTrackingRecHit *hit) const {
0022       std::hash<uint32_t> hasher;
0023       return hasher(hit->geographicalId().rawId());
0024     }
0025   };
0026
0027   using RecHitMap =
0028       cmsutil::SimpleAllocHashMultiMap<const TransientTrackingRecHit *, Trajectory *, HashByDetId, EqualsBySharesInput>;
0029   using TrajMap = cmsutil::UnsortedDumbVectorMap<Trajectory *, int>;
0030
0031   struct Maps {
0032     Maps() : theRecHitMap(128, 256, 1024) {}  // allocate 128 buckets, one row for 256 keys and one row for 512 values
0033     RecHitMap theRecHitMap;
0034     TrajMap theTrajMap;
0035   };
0036
0037   thread_local Maps theMaps;
0038 }  // namespace
0039
0040 using namespace std;
0041
0042 void TrajectoryCleanerBySharedHits::clean(TrajectoryPointerContainer &tc) const {
0043   if (tc.size() <= 1)
0044     return;  // nothing to clean
0045
0046   auto &theRecHitMap = theMaps.theRecHitMap;
0047
0048   theRecHitMap.clear(10 * tc.size());  // set 10*tc.size() active buckets
0049                                        // numbers are not optimized
0050
0051   DEBUG_PRINT(std::cout << "Filling RecHit map" << std::endl);
0052   for (auto const &it : tc) {
0053     DEBUG_PRINT(std::cout << "  Processing trajectory " << it << " (" << it->foundHits() << " valid hits)"
0054                           << std::endl);
0055     auto const &pd = it->measurements();
0056     for (auto const &im : pd) {
0057       auto theRecHit = &(*im.recHit());
0058       if (theRecHit->isValid()) {
0059         DEBUG_PRINT(std::cout << "    Added hit " << theRecHit << " for trajectory " << it << std::endl);
0060         theRecHitMap.insert(theRecHit, it);
0061       }
0062     }
0063   }
0064   DEBUG_PRINT(theRecHitMap.dump());
0065
0066   DEBUG_PRINT(std::cout << "Using RecHit map" << std::endl);
0067   // for each trajectory fill theTrajMap
0068   auto &theTrajMap = theMaps.theTrajMap;
0069   for (auto const &itt : tc) {
0070     if (itt->isValid()) {
0071       DEBUG_PRINT(std::cout << "  Processing trajectory " << itt << " (" << itt->foundHits() << " valid hits)"
0072                             << std::endl);
0073       theTrajMap.clear();
0074       const Trajectory::DataContainer &pd = itt->measurements();
0075       for (auto const &im : pd) {
0076         auto theRecHit = &(*im.recHit());
0077         if (theRecHit->isValid()) {
0078           DEBUG_PRINT(std::cout << "    Searching for overlaps on hit " << theRecHit << " for trajectory " << itt
0079                                 << std::endl);
0080           for (RecHitMap::value_iterator ivec = theRecHitMap.values(theRecHit); ivec.good(); ++ivec) {
0081             if (*ivec != itt) {
0082               if ((*ivec)->isValid()) {
0083                 theTrajMap[*ivec]++;
0084               }
0085             }
0086           }
0087         }
0088       }
0089       //end filling theTrajMap
0090
0091       auto score = [&](Trajectory const &t) -> float {
0092         // possible variant under study
0093         // auto ns = t.foundHits()-t.trailingFoundHits();
0094         //auto penalty =  0.8f*missingHitPenalty_;
0095         // return validHitBonus_*(t.foundHits()-0.2f*t.cccBadHits())  - penalty*t.lostHits() - t.chiSquared();
0096         // classical score
0097         return validHitBonus_ * t.foundHits() - missingHitPenalty_ * t.lostHits() - t.chiSquared();
0098       };
0099
0100       // check for duplicated tracks
0101       if (!theTrajMap.empty() > 0) {
0102         for (auto const &imapp : theTrajMap) {
0103           if (imapp.second > 0) {  // at least 1 hits in common!!!
0104             int innerHit = 0;
0105             if (allowSharedFirstHit) {
0106               const TrajectoryMeasurement &innerMeasure1 =
0107                   (itt->direction() == alongMomentum) ? itt->firstMeasurement() : itt->lastMeasurement();
0108               const TransientTrackingRecHit *h1 = &(*(innerMeasure1).recHit());
0109               const TrajectoryMeasurement &innerMeasure2 = (imapp.first->direction() == alongMomentum)
0110                                                                ? imapp.first->firstMeasurement()
0111                                                                : imapp.first->lastMeasurement();
0112               const TransientTrackingRecHit *h2 = &(*(innerMeasure2).recHit());
0113               if ((h1 == h2) || ((h1->geographicalId() == h2->geographicalId()) &&
0114                                  (h1->hit()->sharesInput(h2->hit(), TrackingRecHit::some)))) {
0115                 innerHit = 1;
0116               }
0117             }
0118             int nhit1 = itt->foundHits();
0119             int nhit2 = imapp.first->foundHits();
0120             if ((imapp.second - innerHit) >= ((min(nhit1, nhit2) - innerHit) * theFraction)) {
0121               Trajectory *badtraj;
0122               auto score1 = score(*itt);
0123               auto score2 = score(*imapp.first);
0124               badtraj = (score1 > score2) ? imapp.first : itt;
0125               badtraj->invalidate();  // invalidate this trajectory
0126             }
0127           }
0128         }
0129       }
0130     }
0131   }
0132 }