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