TrackReco/src/HitPattern.cc

0001 #include "DataFormats/TrackReco/interface/HitPattern.h"
0002 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0003 #include "DataFormats/MuonDetId/interface/DTLayerId.h"
0004 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0005 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0006 #include "DataFormats/MuonDetId/interface/GEMDetId.h"
0007 #include "DataFormats/MuonDetId/interface/ME0DetId.h"
0008 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0009 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0010
0011 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0012
0013 #include "FWCore/Utilities/interface/Likely.h"
0014
0015 #include <bitset>
0016
0017 using namespace reco;
0018
0019 HitPattern::HitPattern()
0020     : hitCount(0), beginTrackHits(0), endTrackHits(0), beginInner(0), endInner(0), beginOuter(0), endOuter(0) {
0021   memset(hitPattern, HitPattern::EMPTY_PATTERN, sizeof(uint16_t) * HitPattern::ARRAY_LENGTH);
0022 }
0023
0024 HitPattern::HitPattern(const HitPattern& other)
0025     : hitCount(other.hitCount),
0026       beginTrackHits(other.beginTrackHits),
0027       endTrackHits(other.endTrackHits),
0028       beginInner(other.beginInner),
0029       endInner(other.endInner),
0030       beginOuter(other.beginOuter),
0031       endOuter(other.endOuter) {
0032   memcpy(this->hitPattern, other.hitPattern, sizeof(uint16_t) * HitPattern::ARRAY_LENGTH);
0033 }
0034
0035 HitPattern::HitPattern(const Run3ScoutingHitPatternPOD& other)
0036     : hitCount(other.hitCount),
0037       beginTrackHits(other.beginTrackHits),
0038       endTrackHits(other.endTrackHits),
0039       beginInner(other.beginInner),
0040       endInner(other.endInner),
0041       beginOuter(other.beginOuter),
0042       endOuter(other.endOuter) {
0043   const unsigned short max_vector_length =
0044       (other.hitPattern.size() > HitPattern::ARRAY_LENGTH) ? HitPattern::ARRAY_LENGTH : other.hitPattern.size();
0045   std::copy(other.hitPattern.begin(), other.hitPattern.begin() + max_vector_length, this->hitPattern);
0046 }
0047
0048 HitPattern::~HitPattern() { ; }
0049
0050 HitPattern& HitPattern::operator=(const HitPattern& other) {
0051   if (this == &other) {
0052     return *this;
0053   }
0054
0055   this->hitCount = other.hitCount;
0056
0057   this->beginTrackHits = other.beginTrackHits;
0058   this->endTrackHits = other.endTrackHits;
0059
0060   this->beginInner = other.beginInner;
0061   this->endInner = other.endInner;
0062
0063   this->beginOuter = other.beginOuter;
0064   this->endOuter = other.endOuter;
0065
0066   memcpy(this->hitPattern, other.hitPattern, sizeof(uint16_t) * HitPattern::ARRAY_LENGTH);
0067
0068   return *this;
0069 }
0070
0071 void HitPattern::clear(void) {
0072   this->hitCount = 0;
0073   this->beginTrackHits = 0;
0074   this->endTrackHits = 0;
0075   this->beginInner = 0;
0076   this->endInner = 0;
0077   this->beginOuter = 0;
0078   this->endOuter = 0;
0079
0080   memset(this->hitPattern, EMPTY_PATTERN, sizeof(uint16_t) * HitPattern::ARRAY_LENGTH);
0081 }
0082
0083 bool HitPattern::appendHit(const TrackingRecHitRef& ref, const TrackerTopology& ttopo) {
0084   return appendHit(*ref, ttopo);
0085 }
0086
0087 uint16_t HitPattern::encode(const TrackingRecHit& hit, const TrackerTopology& ttopo) {
0088   return encode(hit.geographicalId(), hit.getType(), ttopo);
0089 }
0090
0091 namespace {
0092   uint16_t encodeMuonLayer(const DetId& id) {
0093     uint16_t detid = id.det();
0094     uint16_t subdet = id.subdetId();
0095
0096     uint16_t layer = 0x0;
0097     if (detid == DetId::Muon) {
0098       switch (subdet) {
0099         case MuonSubdetId::DT:
0100           layer = ((DTLayerId(id.rawId()).station() - 1) << 2);
0101           layer |= DTLayerId(id.rawId()).superLayer();
0102           break;
0103         case MuonSubdetId::CSC:
0104           layer = ((CSCDetId(id.rawId()).station() - 1) << 2);
0105           layer |= (CSCDetId(id.rawId()).ring() - 1);
0106           break;
0107         case MuonSubdetId::RPC: {
0108           RPCDetId rpcid(id.rawId());
0109           layer = ((rpcid.station() - 1) << 2);
0110           layer |= (rpcid.station() <= 2) ? ((rpcid.layer() - 1) << 1) : 0x0;
0111           layer |= abs(rpcid.region());
0112         } break;
0113         case MuonSubdetId::GEM: {
0114           GEMDetId gemid(id.rawId());
0115           {
0116             uint16_t st = gemid.station();
0117             uint16_t la = gemid.layer();
0118             if (st == 0) {
0119               layer |= 0b1000;
0120               layer |= (la - 1);
0121             } else {
0122               layer |= (st - 1) << 2;
0123               layer |= (la - 1);
0124             }
0125           }
0126         } break;
0127         case MuonSubdetId::ME0: {
0128           ME0DetId me0id(id.rawId());
0129           //layer = ((me0id.roll()-1)<<1) + abs(me0id.layer()-1);
0130           //layer = ((me0id.roll()-1)<<1) + abs(me0id.layer());
0131           //Only layer information that is meaningful is in the roll/etapartition
0132           layer = (me0id.roll());
0133         } break;
0134       }
0135     }
0136     return layer;
0137   }
0138   uint16_t encodeTimingLayer(const DetId& id) {
0139     uint16_t detid = id.det();
0140     uint16_t subdet = id.subdetId();
0141     uint16_t layer = 0x0;
0142     if (detid == DetId::Forward && subdet == FastTime) {
0143       MTDDetId mtdid(id);
0144       switch (mtdid.mtdSubDetector()) {
0145         case MTDDetId::BTL:
0146           layer = BTLDetId(id).modType();
0147           break;
0148         case MTDDetId::ETL:
0149           layer = ETLDetId(id).mtdRR();
0150           break;
0151         default:
0152           throw cms::Exception("HitPattern") << "Invalid MTD Subdetector " << mtdid.mtdSubDetector() << "!";
0153       }
0154     } else {
0155       throw cms::Exception("HitPattern") << "Invalid DetId for FastTime det= " << detid << " subdet= " << subdet << "!";
0156     }
0157     return layer;
0158   }
0159 }  // namespace
0160
0161 uint16_t HitPattern::encode(const DetId& id, TrackingRecHit::Type hitType, const TrackerTopology& ttopo) {
0162   uint16_t detid = id.det();
0163   uint16_t subdet = id.subdetId();
0164
0165   // adding layer/disk/wheel bits
0166   uint16_t layer = 0x0;
0167   if (detid == DetId::Tracker) {
0168     layer = ttopo.layer(id);
0169   } else if (detid == DetId::Muon) {
0170     layer = encodeMuonLayer(id);
0171   } else if (detid == DetId::Forward && subdet == FastTime) {
0172     layer = encodeTimingLayer(id);
0173   }
0174
0175   // adding mono/stereo bit
0176   uint16_t side = 0x0;
0177   if (detid == DetId::Tracker) {
0178     side = isStereo(id, ttopo);
0179   } else if (detid == DetId::Muon || (detid == DetId::Forward && subdet == FastTime)) {
0180     side = 0x0;
0181   }
0182
0183   // juggle the detid around to deal with the fact the bitwidth is larger
0184   // DetId::Muon is 2 and DetId::Forward is 6, must map to 0 and 2 respectively
0185   if (detid == DetId::Tracker) {
0186     detid = TRACKER_HIT;
0187   } else if (detid == DetId::Muon) {
0188     detid = MUON_HIT;  // DetId::Muon is 2 and needs to be reordered to match old encoding where it got masked
0189   } else if (detid == DetId::Forward && subdet == FastTime) {
0190     detid = MTD_HIT;  // since DetId::Forward is some other number, reorder it here
0191   }
0192
0193   return encode(detid, subdet, layer, side, hitType);
0194 }
0195
0196 uint16_t HitPattern::encode(uint16_t det, uint16_t subdet, uint16_t layer, uint16_t side, TrackingRecHit::Type hitType) {
0197   uint16_t pattern = HitPattern::EMPTY_PATTERN;
0198
0199   // adding tracker/muon/mtd detector bits
0200   pattern |= (det & SubDetectorMask) << SubDetectorOffset;
0201
0202   // adding substructure (PXB, PXF, TIB, TID, TOB, TEC, or DT, CSC, RPC,GEM, or BTL, ETL) bits
0203   pattern |= (subdet & SubstrMask) << SubstrOffset;
0204
0205   // adding layer/disk/wheel/ring/modType bits
0206   pattern |= (layer & LayerMask) << LayerOffset;
0207
0208   // adding mono/stereo bit
0209   pattern |= (side & SideMask) << SideOffset;
0210
0211   TrackingRecHit::Type patternHitType =
0212       (hitType == TrackingRecHit::missing_inner || hitType == TrackingRecHit::missing_outer)
0213           ? TrackingRecHit::missing
0214           : ((hitType == TrackingRecHit::inactive_inner || hitType == TrackingRecHit::inactive_outer)
0215                  ? TrackingRecHit::inactive
0216                  : hitType);
0217
0218   pattern |= (patternHitType & HitTypeMask) << HitTypeOffset;
0219
0220   return pattern;
0221 }
0222
0223 bool HitPattern::appendHit(const TrackingRecHit& hit, const TrackerTopology& ttopo) {
0224   return appendHit(hit.geographicalId(), hit.getType(), ttopo);
0225 }
0226
0227 bool HitPattern::appendHit(const DetId& id, TrackingRecHit::Type hitType, const TrackerTopology& ttopo) {
0228   //if HitPattern is full, journey ends no matter what.
0229   if UNLIKELY ((hitCount == HitPattern::MaxHits)) {
0230     return false;
0231   }
0232
0233   uint16_t pattern = HitPattern::encode(id, hitType, ttopo);
0234
0235   return appendHit(pattern, hitType);
0236 }
0237
0238 bool HitPattern::appendHit(const uint16_t pattern, TrackingRecHit::Type hitType) {
0239   //if HitPattern is full, journey ends no matter what.
0240   if UNLIKELY ((hitCount == HitPattern::MaxHits)) {
0241     return false;
0242   }
0243
0244   switch (hitType) {
0245     case TrackingRecHit::valid:
0246     case TrackingRecHit::missing:
0247     case TrackingRecHit::inactive:
0248     case TrackingRecHit::bad:
0249       // hitCount != endT => we are not inserting T type of hits but of T'
0250       // 0 != beginT || 0 != endT => we already have hits of T type
0251       // so we already have hits of T in the vector and we don't want to
0252       // mess them with T' hits.
0253       if UNLIKELY (((hitCount != endTrackHits) && (0 != beginTrackHits || 0 != endTrackHits))) {
0254         cms::Exception("HitPattern") << "TRACK_HITS"
0255                                      << " were stored on this object before hits of some other category were inserted "
0256                                      << "but hits of the same category should be inserted in a row. "
0257                                      << "Please rework the code so it inserts all "
0258                                      << "TRACK_HITS"
0259                                      << " in a row.";
0260         return false;
0261       }
0262       return insertTrackHit(pattern);
0263       break;
0264     case TrackingRecHit::inactive_inner:
0265     case TrackingRecHit::missing_inner:
0266       if UNLIKELY (((hitCount != endInner) && (0 != beginInner || 0 != endInner))) {
0267         cms::Exception("HitPattern") << "MISSING_INNER_HITS"
0268                                      << " were stored on this object before hits of some other category were inserted "
0269                                      << "but hits of the same category should be inserted in a row. "
0270                                      << "Please rework the code so it inserts all "
0271                                      << "MISSING_INNER_HITS"
0272                                      << " in a row.";
0273         return false;
0274       }
0275       return insertExpectedInnerHit(pattern);
0276       break;
0277     case TrackingRecHit::inactive_outer:
0278     case TrackingRecHit::missing_outer:
0279       if UNLIKELY (((hitCount != endOuter) && (0 != beginOuter || 0 != endOuter))) {
0280         cms::Exception("HitPattern") << "MISSING_OUTER_HITS"
0281                                      << " were stored on this object before hits of some other category were inserted "
0282                                      << "but hits of the same category should be inserted in a row. "
0283                                      << "Please rework the code so it inserts all "
0284                                      << "MISSING_OUTER_HITS"
0285                                      << " in a row.";
0286         return false;
0287       }
0288       return insertExpectedOuterHit(pattern);
0289       break;
0290   }
0291
0292   return false;
0293 }
0294
0295 bool HitPattern::appendTrackerHit(uint16_t subdet, uint16_t layer, uint16_t stereo, TrackingRecHit::Type hitType) {
0296   return appendHit(encode(TRACKER_HIT, subdet, layer, stereo, hitType), hitType);
0297 }
0298
0299 bool HitPattern::appendMuonHit(const DetId& id, TrackingRecHit::Type hitType) {
0300   //if HitPattern is full, journey ends no matter what.
0301   if UNLIKELY ((hitCount == HitPattern::MaxHits)) {
0302     return false;
0303   }
0304
0305   if UNLIKELY (id.det() != DetId::Muon) {
0306     throw cms::Exception("HitPattern")
0307         << "Got DetId from det " << id.det()
0308         << " that is not Muon in appendMuonHit(), which should only be used for muon hits in the HitPattern IO rule";
0309   }
0310
0311   uint16_t subdet = id.subdetId();
0312   return appendHit(encode(MUON_HIT, subdet, encodeMuonLayer(id), 0, hitType), hitType);
0313 }
0314
0315 uint16_t HitPattern::getHitPatternByAbsoluteIndex(int position) const {
0316   if UNLIKELY ((position < 0 || position >= hitCount)) {
0317     return HitPattern::EMPTY_PATTERN;
0318   }
0319   /*
0320     Note: you are not taking a consecutive sequence of HIT_LENGTH bits starting from position * HIT_LENGTH
0321      as the bit order in the words are reversed.
0322      e.g. if position = 0 you take the lowest 12 bits of the first word.
0323
0324      I hope this can clarify what is the memory layout of such thing
0325
0326     straight 01234567890123456789012345678901 | 23456789012345678901234567890123 | 4567
0327     (global) 0         1         2         3  | 3       4         5         6    | 6
0328     words    [--------------0---------------] | [--------------1---------------] | [---
0329     word     01234567890123456789012345678901 | 01234567890123456789012345678901 | 0123
0330     (str)   0         1         2         3  | 0         1         2         3  | 0
0331           [--------------0---------------] | [--------------1---------------] | [---
0332     word     10987654321098765432109876543210 | 10987654321098765432109876543210 | 1098
0333     (rev)   32         21        10        0 | 32         21        10        0 | 32
0334     reverse  10987654321098765432109876543210 | 32109876543210987654321098765432 | 5432
0335              32         21        10        0 | 6  65        54        43      3   9
0336
0337      ugly enough, but it's not my fault, I was not even in CMS at that time   [gpetrucc]
0338     */
0339
0340   uint16_t bitEndOffset = (position + 1) * HIT_LENGTH;
0341   uint8_t secondWord = (bitEndOffset >> 4);
0342   uint8_t secondWordBits = bitEndOffset & (16 - 1);  // that is, bitEndOffset % 16
0343   if (secondWordBits >= HIT_LENGTH) {                // full block is in this word
0344     uint8_t lowBitsToTrash = secondWordBits - HIT_LENGTH;
0345     uint16_t myResult = (hitPattern[secondWord] >> lowBitsToTrash) & ((1 << HIT_LENGTH) - 1);
0346     return myResult;
0347   } else {
0348     uint8_t firstWordBits = HIT_LENGTH - secondWordBits;
0349     uint16_t firstWordBlock = hitPattern[secondWord - 1] >> (16 - firstWordBits);
0350     if (secondWordBits == 0)
0351       return firstWordBlock;
0352     uint16_t secondWordBlock = hitPattern[secondWord] & ((1 << secondWordBits) - 1);
0353     uint16_t myResult = firstWordBlock + (secondWordBlock << firstWordBits);
0354     return myResult;
0355   }
0356 }
0357
0358 bool HitPattern::hasValidHitInPixelLayer(enum PixelSubdetector::SubDetector det, uint16_t layer) const {
0359   for (int i = beginTrackHits; i < endTrackHits; ++i) {
0360     uint16_t pattern = getHitPatternByAbsoluteIndex(i);
0361     bool pixelHitFilter = ((det == 1 && pixelBarrelHitFilter(pattern)) || (det == 2 && pixelEndcapHitFilter(pattern)));
0362     if (pixelHitFilter && (getLayer(pattern) == layer) && validHitFilter(pattern)) {
0363       return true;
0364     }
0365   }
0366   return false;
0367 }
0368
0369 int HitPattern::numberOfValidStripLayersWithMonoAndStereo(uint16_t stripdet, uint16_t layer) const {
0370   bool hasMono[SubstrMask + 1][LayerMask + 1];
0371   bool hasStereo[SubstrMask + 1][LayerMask + 1];
0372   memset(hasMono, 0, sizeof(hasMono));
0373   memset(hasStereo, 0, sizeof(hasStereo));
0374
0375   // mark which layers have mono/stereo hits
0376   for (int i = beginTrackHits; i < endTrackHits; ++i) {
0377     uint16_t pattern = getHitPatternByAbsoluteIndex(i);
0378     uint16_t subStructure = getSubStructure(pattern);
0379
0380     if (validHitFilter(pattern) && stripHitFilter(pattern)) {
0381       if (stripdet != 0 && subStructure != stripdet) {
0382         continue;
0383       }
0384
0385       if (layer != 0 && getSubSubStructure(pattern) != layer) {
0386         continue;
0387       }
0388
0389       switch (getSide(pattern)) {
0390         case 0:  // mono
0391           hasMono[subStructure][getLayer(pattern)] = true;
0392           break;
0393         case 1:  // stereo
0394           hasStereo[subStructure][getLayer(pattern)] = true;
0395           break;
0396         default:;
0397           break;
0398       }
0399     }
0400   }
0401
0402   // count how many layers have mono and stereo hits
0403   int count = 0;
0404   for (int i = 0; i < SubstrMask + 1; ++i) {
0405     for (int j = 0; j < LayerMask + 1; ++j) {
0406       if (hasMono[i][j] && hasStereo[i][j]) {
0407         count++;
0408       }
0409     }
0410   }
0411   return count;
0412 }
0413
0414 int HitPattern::numberOfValidStripLayersWithMonoAndStereo() const {
0415   auto category = TRACK_HITS;
0416   std::bitset<128> side[2];
0417   std::pair<uint8_t, uint8_t> range = getCategoryIndexRange(category);
0418   for (int i = range.first; i < range.second; ++i) {
0419     auto pattern = getHitPatternByAbsoluteIndex(i);
0420     if (pattern > maxTrackerWord)
0421       continue;
0422     if (pattern < minStripWord)
0423       continue;
0424     uint16_t hitType = (pattern >> HitTypeOffset) & HitTypeMask;
0425     if (hitType != HIT_TYPE::VALID)
0426       continue;
0427     auto apattern = (pattern - minTrackerWord) >> LayerOffset;
0428     // assert(apattern<128);
0429     side[getSide(pattern)].set(apattern);
0430   }
0431   // assert(numberOfValidStripLayersWithMonoAndStereo(0, 0)==int((side[0]&side[1]).count()));
0432   return (side[0] & side[1]).count();
0433 }
0434
0435 int HitPattern::numberOfValidTOBLayersWithMonoAndStereo(uint32_t layer) const {
0436   return numberOfValidStripLayersWithMonoAndStereo(StripSubdetector::TOB, layer);
0437 }
0438
0439 int HitPattern::numberOfValidTIBLayersWithMonoAndStereo(uint32_t layer) const {
0440   return numberOfValidStripLayersWithMonoAndStereo(StripSubdetector::TIB, layer);
0441 }
0442
0443 int HitPattern::numberOfValidTIDLayersWithMonoAndStereo(uint32_t layer) const {
0444   return numberOfValidStripLayersWithMonoAndStereo(StripSubdetector::TID, layer);
0445 }
0446
0447 int HitPattern::numberOfValidTECLayersWithMonoAndStereo(uint32_t layer) const {
0448   return numberOfValidStripLayersWithMonoAndStereo(StripSubdetector::TEC, layer);
0449 }
0450
0451 uint32_t HitPattern::getTrackerLayerCase(HitCategory category, uint16_t substr, uint16_t layer) const {
0452   uint16_t tk_substr_layer =
0453       (0x1 << SubDetectorOffset) + ((substr & SubstrMask) << SubstrOffset) + ((layer & LayerMask) << LayerOffset);
0454
0455   uint16_t mask = (SubDetectorMask << SubDetectorOffset) + (SubstrMask << SubstrOffset) + (LayerMask << LayerOffset);
0456
0457   // layer case 0: valid + (missing, off, bad) ==> with measurement
0458   // layer case 1: missing + (off, bad) ==> without measurement
0459   // layer case 2: off, bad ==> totally off or bad, cannot say much
0460   // layer case NULL_RETURN: track outside acceptance or in gap ==> null
0461   uint32_t layerCase = NULL_RETURN;
0462   std::pair<uint8_t, uint8_t> range = getCategoryIndexRange(category);
0463   for (int i = range.first; i < range.second; ++i) {
0464     uint16_t pattern = getHitPatternByAbsoluteIndex(i);
0465     if ((pattern & mask) == tk_substr_layer) {
0466       uint16_t hitType = (pattern >> HitTypeOffset) & HitTypeMask;
0467       if (hitType < layerCase) {
0468         // BAD and INACTIVE as the same type (as INACTIVE)
0469         layerCase = (hitType == HIT_TYPE::BAD ? (uint32_t)HIT_TYPE::INACTIVE : hitType);
0470         if (layerCase == HIT_TYPE::VALID) {
0471           break;
0472         }
0473       }
0474     }
0475   }
0476   return layerCase;
0477 }
0478
0479 uint16_t HitPattern::getTrackerMonoStereo(HitCategory category, uint16_t substr, uint16_t layer) const {
0480   uint16_t tk_substr_layer =
0481       (0x1 << SubDetectorOffset) + ((substr & SubstrMask) << SubstrOffset) + ((layer & LayerMask) << LayerOffset);
0482   uint16_t mask = (SubDetectorMask << SubDetectorOffset) + (SubstrMask << SubstrOffset) + (LayerMask << LayerOffset);
0483
0484   //             0: neither a valid mono nor a valid stereo hit
0485   //          MONO: valid mono hit
0486   //        STEREO: valid stereo hit
0487   // MONO | STEREO: both
0488   uint16_t monoStereo = 0x0;
0489   std::pair<uint8_t, uint8_t> range = getCategoryIndexRange(category);
0490   for (int i = range.first; i < range.second; ++i) {
0491     uint16_t pattern = getHitPatternByAbsoluteIndex(i);
0492     if ((pattern & mask) == tk_substr_layer) {
0493       uint16_t hitType = (pattern >> HitTypeOffset) & HitTypeMask;
0494       if (hitType == HIT_TYPE::VALID) {
0495         switch (getSide(pattern)) {
0496           case 0:  // mono
0497             monoStereo |= MONO;
0498             break;
0499           case 1:  // stereo
0500             monoStereo |= STEREO;
0501             break;
0502         }
0503       }
0504
0505       if (monoStereo == (MONO | STEREO)) {
0506         break;
0507       }
0508     }
0509   }
0510   return monoStereo;
0511 }
0512
0513 int HitPattern::pixelLayersWithMeasurement() const {
0514   auto category = TRACK_HITS;
0515   std::bitset<128> layerOk;
0516   std::pair<uint8_t, uint8_t> range = getCategoryIndexRange(category);
0517   for (int i = range.first; i < range.second; ++i) {
0518     auto pattern = getHitPatternByAbsoluteIndex(i);
0519     if UNLIKELY (!trackerHitFilter(pattern))
0520       continue;
0521     if (pattern > minStripWord)
0522       continue;
0523     uint16_t hitType = (pattern >> HitTypeOffset) & HitTypeMask;
0524     if (hitType != HIT_TYPE::VALID)
0525       continue;
0526     pattern = (pattern - minTrackerWord) >> LayerOffset;
0527     // assert(pattern<128);
0528     layerOk.set(pattern);
0529   }
0530   // assert(pixelLayersWithMeasurementOld()==int(layerOk.count()));
0531   return layerOk.count();
0532 }
0533
0534 int HitPattern::trackerLayersWithMeasurement() const {
0535   auto category = TRACK_HITS;
0536   std::bitset<128> layerOk;
0537   std::pair<uint8_t, uint8_t> range = getCategoryIndexRange(category);
0538   for (int i = range.first; i < range.second; ++i) {
0539     auto pattern = getHitPatternByAbsoluteIndex(i);
0540     if UNLIKELY (!trackerHitFilter(pattern))
0541       continue;
0542     uint16_t hitType = (pattern >> HitTypeOffset) & HitTypeMask;
0543     if (hitType != HIT_TYPE::VALID)
0544       continue;
0545     pattern = (pattern - minTrackerWord) >> LayerOffset;
0546     // assert(pattern<128);
0547     layerOk.set(pattern);
0548   }
0549   // assert(trackerLayersWithMeasurementOld()==int(layerOk.count()));
0550   return layerOk.count();
0551 }
0552
0553 int HitPattern::trackerLayersWithoutMeasurement(HitCategory category) const {
0554   std::bitset<128> layerOk;
0555   std::bitset<128> layerMissed;
0556   std::pair<uint8_t, uint8_t> range = getCategoryIndexRange(category);
0557   for (int i = range.first; i < range.second; ++i) {
0558     auto pattern = getHitPatternByAbsoluteIndex(i);
0559     if UNLIKELY (!trackerHitFilter(pattern))
0560       continue;
0561     uint16_t hitType = (pattern >> HitTypeOffset) & HitTypeMask;
0562     pattern = (pattern - minTrackerWord) >> LayerOffset;
0563     // assert(pattern<128);
0564     if (hitType == HIT_TYPE::VALID)
0565       layerOk.set(pattern);
0566     if (hitType == HIT_TYPE::MISSING)
0567       layerMissed.set(pattern);
0568   }
0569   layerMissed &= ~layerOk;
0570
0571   // assert(trackerLayersWithoutMeasurementOld(category)==int(layerMissed.count()));
0572
0573   return layerMissed.count();
0574 }
0575
0576 int HitPattern::pixelBarrelLayersWithMeasurement() const {
0577   int count = 0;
0578   uint16_t NPixBarrel = 4;
0579   for (uint16_t layer = 1; layer <= NPixBarrel; layer++) {
0580     if (getTrackerLayerCase(TRACK_HITS, PixelSubdetector::PixelBarrel, layer) == HIT_TYPE::VALID) {
0581       count++;
0582     }
0583   }
0584   return count;
0585 }
0586
0587 int HitPattern::pixelEndcapLayersWithMeasurement() const {
0588   int count = 0;
0589   uint16_t NPixForward = 3;
0590   for (uint16_t layer = 1; layer <= NPixForward; layer++) {
0591     if (getTrackerLayerCase(TRACK_HITS, PixelSubdetector::PixelEndcap, layer) == HIT_TYPE::VALID) {
0592       count++;
0593     }
0594   }
0595   return count;
0596 }
0597
0598 int HitPattern::stripTIBLayersWithMeasurement() const {
0599   int count = 0;
0600   for (uint16_t layer = 1; layer <= 4; layer++) {
0601     if (getTrackerLayerCase(TRACK_HITS, StripSubdetector::TIB, layer) == HIT_TYPE::VALID) {
0602       count++;
0603     }
0604   }
0605   return count;
0606 }
0607
0608 int HitPattern::stripTIDLayersWithMeasurement() const {
0609   int count = 0;
0610   for (uint16_t layer = 1; layer <= 3; layer++) {
0611     if (getTrackerLayerCase(TRACK_HITS, StripSubdetector::TID, layer) == HIT_TYPE::VALID) {
0612       count++;
0613     }
0614   }
0615   return count;
0616 }
0617
0618 int HitPattern::stripTOBLayersWithMeasurement() const {
0619   int count = 0;
0620   for (uint16_t layer = 1; layer <= 6; layer++) {
0621     if (getTrackerLayerCase(TRACK_HITS, StripSubdetector::TOB, layer) == HIT_TYPE::VALID) {
0622       count++;
0623     }
0624   }
0625   return count;
0626 }
0627
0628 int HitPattern::stripTECLayersWithMeasurement() const {
0629   int count = 0;
0630   for (uint16_t layer = 1; layer <= 9; layer++) {
0631     if (getTrackerLayerCase(TRACK_HITS, StripSubdetector::TEC, layer) == HIT_TYPE::VALID) {
0632       count++;
0633     }
0634   }
0635   return count;
0636 }
0637
0638 int HitPattern::pixelBarrelLayersWithoutMeasurement(HitCategory category) const {
0639   int count = 0;
0640   uint16_t NPixBarrel = 4;
0641   for (uint16_t layer = 1; layer <= NPixBarrel; layer++) {
0642     if (getTrackerLayerCase(category, PixelSubdetector::PixelBarrel, layer) == HIT_TYPE::MISSING) {
0643       count++;
0644     }
0645   }
0646   return count;
0647 }
0648
0649 int HitPattern::pixelEndcapLayersWithoutMeasurement(HitCategory category) const {
0650   int count = 0;
0651   uint16_t NPixForward = 3;
0652   for (uint16_t layer = 1; layer <= NPixForward; layer++) {
0653     if (getTrackerLayerCase(category, PixelSubdetector::PixelEndcap, layer) == HIT_TYPE::MISSING) {
0654       count++;
0655     }
0656   }
0657   return count;
0658 }
0659
0660 int HitPattern::stripTIBLayersWithoutMeasurement(HitCategory category) const {
0661   int count = 0;
0662   for (uint16_t layer = 1; layer <= 4; layer++) {
0663     if (getTrackerLayerCase(category, StripSubdetector::TIB, layer) == HIT_TYPE::MISSING) {
0664       count++;
0665     }
0666   }
0667   return count;
0668 }
0669
0670 int HitPattern::stripTIDLayersWithoutMeasurement(HitCategory category) const {
0671   int count = 0;
0672   for (uint16_t layer = 1; layer <= 3; layer++) {
0673     if (getTrackerLayerCase(category, StripSubdetector::TID, layer) == HIT_TYPE::MISSING) {
0674       count++;
0675     }
0676   }
0677   return count;
0678 }
0679
0680 int HitPattern::stripTOBLayersWithoutMeasurement(HitCategory category) const {
0681   int count = 0;
0682   for (uint16_t layer = 1; layer <= 6; layer++) {
0683     if (getTrackerLayerCase(category, StripSubdetector::TOB, layer) == HIT_TYPE::MISSING) {
0684       count++;
0685     }
0686   }
0687   return count;
0688 }
0689
0690 int HitPattern::stripTECLayersWithoutMeasurement(HitCategory category) const {
0691   int count = 0;
0692   for (uint16_t layer = 1; layer <= 9; layer++) {
0693     if (getTrackerLayerCase(category, StripSubdetector::TEC, layer) == HIT_TYPE::MISSING) {
0694       count++;
0695     }
0696   }
0697   return count;
0698 }
0699
0700 int HitPattern::pixelBarrelLayersTotallyOffOrBad(HitCategory category) const {
0701   int count = 0;
0702   uint16_t NPixBarrel = 4;
0703   for (uint16_t layer = 1; layer <= NPixBarrel; layer++) {
0704     if (getTrackerLayerCase(category, PixelSubdetector::PixelBarrel, layer) == HIT_TYPE::INACTIVE) {
0705       count++;
0706     }
0707   }
0708   return count;
0709 }
0710
0711 int HitPattern::pixelEndcapLayersTotallyOffOrBad(HitCategory category) const {
0712   int count = 0;
0713   uint16_t NPixForward = 3;
0714   for (uint16_t layer = 1; layer <= NPixForward; layer++) {
0715     if (getTrackerLayerCase(category, PixelSubdetector::PixelEndcap, layer) == HIT_TYPE::INACTIVE) {
0716       count++;
0717     }
0718   }
0719   return count;
0720 }
0721
0722 int HitPattern::stripTIBLayersTotallyOffOrBad(HitCategory category) const {
0723   int count = 0;
0724   for (uint16_t layer = 1; layer <= 4; layer++) {
0725     if (getTrackerLayerCase(category, StripSubdetector::TIB, layer) == HIT_TYPE::INACTIVE) {
0726       count++;
0727     }
0728   }
0729   return count;
0730 }
0731
0732 int HitPattern::stripTIDLayersTotallyOffOrBad(HitCategory category) const {
0733   int count = 0;
0734   for (uint16_t layer = 1; layer <= 3; layer++) {
0735     if (getTrackerLayerCase(category, StripSubdetector::TID, layer) == HIT_TYPE::INACTIVE) {
0736       count++;
0737     }
0738   }
0739   return count;
0740 }
0741
0742 int HitPattern::stripTOBLayersTotallyOffOrBad(HitCategory category) const {
0743   int count = 0;
0744   for (uint16_t layer = 1; layer <= 6; layer++) {
0745     if (getTrackerLayerCase(category, StripSubdetector::TOB, layer) == HIT_TYPE::INACTIVE) {
0746       count++;
0747     }
0748   }
0749   return count;
0750 }
0751
0752 int HitPattern::stripTECLayersTotallyOffOrBad(HitCategory category) const {
0753   int count = 0;
0754   for (uint16_t layer = 1; layer <= 9; layer++) {
0755     if (getTrackerLayerCase(category, StripSubdetector::TEC, layer) == HIT_TYPE::INACTIVE) {
0756       count++;
0757     }
0758   }
0759   return count;
0760 }
0761
0762 int HitPattern::pixelBarrelLayersNull() const {
0763   int count = 0;
0764   uint16_t NPixBarrel = 4;
0765   for (uint16_t layer = 1; layer <= NPixBarrel; layer++) {
0766     if (getTrackerLayerCase(TRACK_HITS, PixelSubdetector::PixelBarrel, layer) == NULL_RETURN) {
0767       count++;
0768     }
0769   }
0770   return count;
0771 }
0772
0773 int HitPattern::pixelEndcapLayersNull() const {
0774   int count = 0;
0775   uint16_t NPixForward = 3;
0776   for (uint16_t layer = 1; layer <= NPixForward; layer++) {
0777     if (getTrackerLayerCase(TRACK_HITS, PixelSubdetector::PixelEndcap, layer) == NULL_RETURN) {
0778       count++;
0779     }
0780   }
0781   return count;
0782 }
0783
0784 int HitPattern::stripTIBLayersNull() const {
0785   int count = 0;
0786   for (uint16_t layer = 1; layer <= 4; layer++) {
0787     if (getTrackerLayerCase(TRACK_HITS, StripSubdetector::TIB, layer) == NULL_RETURN) {
0788       count++;
0789     }
0790   }
0791   return count;
0792 }
0793
0794 int HitPattern::stripTIDLayersNull() const {
0795   int count = 0;
0796   for (uint16_t layer = 1; layer <= 3; layer++) {
0797     if (getTrackerLayerCase(TRACK_HITS, StripSubdetector::TID, layer) == NULL_RETURN) {
0798       count++;
0799     }
0800   }
0801   return count;
0802 }
0803
0804 int HitPattern::stripTOBLayersNull() const {
0805   int count = 0;
0806   for (uint16_t layer = 1; layer <= 6; layer++) {
0807     if (getTrackerLayerCase(TRACK_HITS, StripSubdetector::TOB, layer) == NULL_RETURN) {
0808       count++;
0809     }
0810   }
0811   return count;
0812 }
0813
0814 int HitPattern::stripTECLayersNull() const {
0815   int count = 0;
0816   for (uint16_t layer = 1; layer <= 9; layer++) {
0817     if (getTrackerLayerCase(TRACK_HITS, StripSubdetector::TEC, layer) == NULL_RETURN) {
0818       count++;
0819     }
0820   }
0821   return count;
0822 }
0823
0824 void HitPattern::printHitPattern(HitCategory category, int position, std::ostream& stream) const {
0825   uint16_t pattern = getHitPattern(category, position);
0826   stream << "\t";
0827   if (muonHitFilter(pattern)) {
0828     stream << "muon";
0829   } else if (trackerHitFilter(pattern)) {
0830     stream << "tracker";
0831   } else if (timingHitFilter(pattern)) {
0832     stream << "timing";
0833   }
0834
0835   stream << "\tsubstructure " << getSubStructure(pattern);
0836   if (muonHitFilter(pattern)) {
0837     stream << "\tstation " << getMuonStation(pattern);
0838     if (muonDTHitFilter(pattern)) {
0839       stream << "\tdt superlayer " << getDTSuperLayer(pattern);
0840     } else if (muonCSCHitFilter(pattern)) {
0841       stream << "\tcsc ring " << getCSCRing(pattern);
0842     } else if (muonRPCHitFilter(pattern)) {
0843       stream << "\trpc " << (getRPCregion(pattern) ? "endcaps" : "barrel") << ", layer " << getRPCLayer(pattern);
0844     } else if (muonGEMHitFilter(pattern)) {
0845       stream << "\tgem "
0846              << " station " << getGEMStation(pattern) << ", layer" << getGEMLayer(pattern);
0847     } else if (muonME0HitFilter(pattern)) {
0848       stream << "\tme0 ";
0849     } else {
0850       stream << "(UNKNOWN Muon SubStructure!) \tsubsubstructure " << getSubStructure(pattern);
0851     }
0852   } else if (timingHitFilter(pattern)) {
0853     stream << "\tdetector " << getSubStructure(pattern);
0854   } else {
0855     stream << "\tlayer " << getLayer(pattern);
0856   }
0857   stream << "\thit type " << getHitType(pattern);
0858   stream << std::endl;
0859 }
0860
0861 void HitPattern::print(HitCategory category, std::ostream& stream) const {
0862   stream << "HitPattern" << std::endl;
0863   for (int i = 0; i < numberOfAllHits(category); ++i) {
0864     printHitPattern(category, i, stream);
0865   }
0866   std::ios_base::fmtflags flags = stream.flags();
0867   stream.setf(std::ios_base::hex, std::ios_base::basefield);
0868   stream.setf(std::ios_base::showbase);
0869
0870   for (int i = 0; i < this->numberOfAllHits(category); ++i) {
0871     stream << getHitPattern(category, i) << std::endl;
0872   }
0873
0874   stream.flags(flags);
0875 }
0876
0877 uint16_t HitPattern::isStereo(DetId i, const TrackerTopology& ttopo) {
0878   if (i.det() != DetId::Tracker) {
0879     return 0;
0880   }
0881
0882   switch (i.subdetId()) {
0883     case PixelSubdetector::PixelBarrel:
0884     case PixelSubdetector::PixelEndcap:
0885       return 0;
0886     case StripSubdetector::TIB:
0887       return ttopo.tibIsStereo(i);
0888     case StripSubdetector::TID:
0889       return ttopo.tidIsStereo(i);
0890     case StripSubdetector::TOB:
0891       return ttopo.tobIsStereo(i);
0892     case StripSubdetector::TEC:
0893       return ttopo.tecIsStereo(i);
0894     default:
0895       return 0;
0896   }
0897 }
0898
0899 int HitPattern::muonStations(int subdet, int hitType) const {
0900   int stations[5] = {0, 0, 0, 0, 0};
0901   for (int i = beginTrackHits; i < endTrackHits; ++i) {
0902     uint16_t pattern = getHitPatternByAbsoluteIndex(i);
0903     if (muonHitFilter(pattern) && (subdet == 0 || int(getSubStructure(pattern)) == subdet) &&
0904         (hitType == -1 || int(getHitType(pattern)) == hitType)) {
0905       stations[getMuonStation(pattern)] = 1;
0906     }
0907   }
0908
0909   return stations[0] + stations[1] + stations[2] + stations[3] + stations[4];
0910 }
0911
0912 int HitPattern::innermostMuonStationWithHits(int hitType) const {
0913   int ret = 0;
0914   for (int i = beginTrackHits; i < endTrackHits; ++i) {
0915     uint16_t pattern = getHitPatternByAbsoluteIndex(i);
0916     if (muonHitFilter(pattern) && (hitType == -1 || int(getHitType(pattern)) == hitType)) {
0917       int stat = getMuonStation(pattern);
0918       if (ret == 0 || stat < ret) {
0919         ret = stat;
0920       }
0921     }
0922   }
0923
0924   return ret;
0925 }
0926
0927 int HitPattern::outermostMuonStationWithHits(int hitType) const {
0928   int ret = 0;
0929   for (int i = beginTrackHits; i < endTrackHits; ++i) {
0930     uint16_t pattern = getHitPatternByAbsoluteIndex(i);
0931     if (muonHitFilter(pattern) && (hitType == -1 || int(getHitType(pattern)) == hitType)) {
0932       int stat = getMuonStation(pattern);
0933       if (ret == 0 || stat > ret) {
0934         ret = stat;
0935       }
0936     }
0937   }
0938   return ret;
0939 }
0940
0941 int HitPattern::numberOfDTStationsWithRPhiView() const {
0942   int stations[4] = {0, 0, 0, 0};
0943   for (int i = beginTrackHits; i < endTrackHits; ++i) {
0944     uint16_t pattern = getHitPatternByAbsoluteIndex(i);
0945
0946     if (muonDTHitFilter(pattern) && validHitFilter(pattern) && getDTSuperLayer(pattern) != 2) {
0947       stations[getMuonStation(pattern) - 1] = 1;
0948     }
0949   }
0950   return stations[0] + stations[1] + stations[2] + stations[3];
0951 }
0952
0953 int HitPattern::numberOfDTStationsWithRZView() const {
0954   int stations[4] = {0, 0, 0, 0};
0955   for (int i = beginTrackHits; i < endTrackHits; ++i) {
0956     uint16_t pattern = getHitPatternByAbsoluteIndex(i);
0957     if (muonDTHitFilter(pattern) && validHitFilter(pattern) && getDTSuperLayer(pattern) == 2) {
0958       stations[getMuonStation(pattern) - 1] = 1;
0959     }
0960   }
0961   return stations[0] + stations[1] + stations[2] + stations[3];
0962 }
0963
0964 int HitPattern::numberOfDTStationsWithBothViews() const {
0965   int stations[4][2] = {{0, 0}, {0, 0}, {0, 0}, {0, 0}};
0966   for (int i = beginTrackHits; i < endTrackHits; ++i) {
0967     uint16_t pattern = getHitPatternByAbsoluteIndex(i);
0968     if (muonDTHitFilter(pattern) && validHitFilter(pattern)) {
0969       stations[getMuonStation(pattern) - 1][getDTSuperLayer(pattern) == 2] = 1;
0970     }
0971   }
0972
0973   return stations[0][0] * stations[0][1] + stations[1][0] * stations[1][1] + stations[2][0] * stations[2][1] +
0974          stations[3][0] * stations[3][1];
0975 }
0976
0977 void HitPattern::insertHit(const uint16_t pattern) {
0978   int offset = hitCount * HIT_LENGTH;
0979   for (int i = 0; i < HIT_LENGTH; i++) {
0980     int pos = offset + i;
0981     uint16_t bit = (pattern >> i) & 0x1;
0982     //equivalent to hitPattern[pos / 16] += bit << ((offset + i) % 16);
0983     hitPattern[pos >> 4] += bit << ((offset + i) & (16 - 1));
0984   }
0985   hitCount++;
0986 }
0987
0988 bool HitPattern::insertTrackHit(const uint16_t pattern) {
0989   // if begin is 0, this is the first hit of this type being inserted, so
0990   // we need to update begin so it points to the correct index, the first
0991   // empty index.
0992   // unlikely, because it will happen only when inserting
0993   // the first hit of this type
0994   if UNLIKELY ((0 == beginTrackHits && 0 == endTrackHits)) {
0995     beginTrackHits = hitCount;
0996     // before the first hit of this type is inserted, there are no hits
0997     endTrackHits = beginTrackHits;
0998   }
0999
1000   insertHit(pattern);
1001   endTrackHits++;
1002
1003   return true;
1004 }
1005
1006 bool HitPattern::insertExpectedInnerHit(const uint16_t pattern) {
1007   if UNLIKELY ((0 == beginInner && 0 == endInner)) {
1008     beginInner = hitCount;
1009     endInner = beginInner;
1010   }
1011
1012   insertHit(pattern);
1013   endInner++;
1014
1015   return true;
1016 }
1017
1018 bool HitPattern::insertExpectedOuterHit(const uint16_t pattern) {
1019   if UNLIKELY ((0 == beginOuter && 0 == endOuter)) {
1020     beginOuter = hitCount;
1021     endOuter = beginOuter;
1022   }
1023
1024   insertHit(pattern);
1025   endOuter++;
1026
1027   return true;
1028 }
1029
1030 Run3ScoutingHitPatternPOD HitPattern::run3ScoutingHitPatternPOD() const {
1031   Run3ScoutingHitPatternPOD result{
1032       .hitCount = hitCount,
1033       .beginTrackHits = beginTrackHits,
1034       .endTrackHits = endTrackHits,
1035       .beginInner = beginInner,
1036       .endInner = endInner,
1037       .beginOuter = beginOuter,
1038       .endOuter = endOuter,
1039       .hitPattern = std::vector<uint16_t>(hitPattern, hitPattern + HitPattern::ARRAY_LENGTH)};
1040   return result;
1041 }