File indexing completed on 2024-04-06 12:31:40
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
0008 #define DEBUG_PRINT(X)
0009
0010 namespace {
0011
0012
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
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) {}
0033 RecHitMap theRecHitMap;
0034 TrajMap theTrajMap;
0035 };
0036
0037 thread_local Maps theMaps;
0038 }
0039
0040 using namespace std;
0041
0042 void TrajectoryCleanerBySharedHits::clean(TrajectoryPointerContainer &tc) const {
0043 if (tc.size() <= 1)
0044 return;
0045
0046 auto &theRecHitMap = theMaps.theRecHitMap;
0047
0048 theRecHitMap.clear(10 * tc.size());
0049
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
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
0090
0091 auto score = [&](Trajectory const &t) -> float {
0092
0093
0094
0095
0096
0097 return validHitBonus_ * t.foundHits() - missingHitPenalty_ * t.lostHits() - t.chiSquared();
0098 };
0099
0100
0101 if (!theTrajMap.empty() > 0) {
0102 for (auto const &imapp : theTrajMap) {
0103 if (imapp.second > 0) {
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();
0126 }
0127 }
0128 }
0129 }
0130 }
0131 }
0132 }