File indexing completed on 2025-02-05 23:51:52
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include <algorithm>
0010 #include <memory>
0011
0012 #include <vector>
0013
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015
0016 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0017 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0018
0019 #include "DataFormats/DetId/interface/DetId.h"
0020 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0021 #include "DataFormats/MuonDetId/interface/DTLayerId.h"
0022 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0023 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0024 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0025 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0026 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0027 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0028
0029 #include "DataFormats/TrackReco/interface/Track.h"
0030 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0031 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0032 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0033
0034 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0035
0036 #include "SimTracker/TrackHistory/interface/TrackQuality.h"
0037
0038
0039
0040 namespace {
0041 const uint32_t NonMatchedTrackId = (uint32_t)-1;
0042
0043 struct MatchedHit {
0044 DetId detId;
0045 uint32_t simTrackId;
0046 EncodedEventId collision;
0047 int recHitId;
0048 TrackQuality::Layer::State state;
0049
0050 bool operator<(const MatchedHit &other) const {
0051 if (detId < other.detId)
0052 return true;
0053 else if (detId > other.detId)
0054 return false;
0055 else if (collision < other.collision)
0056 return true;
0057 else if (other.collision < collision)
0058 return false;
0059 else
0060 return simTrackId < other.simTrackId;
0061 }
0062
0063 bool operator==(const MatchedHit &other) const {
0064 return detId == other.detId && collision == other.collision && simTrackId == other.simTrackId;
0065 }
0066 };
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077 }
0078
0079 typedef std::pair<TrackQuality::Layer::SubDet, short int> DetLayer;
0080
0081
0082 static const int statePriorities[] = {
0083 3,
0084 5,
0085 0,
0086 7,
0087 2,
0088 4,
0089 6,
0090 1};
0091
0092 DetLayer getDetLayer(DetId detId, const TrackerTopology *tTopo) {
0093 TrackQuality::Layer::SubDet det = TrackQuality::Layer::Invalid;
0094 short int layer = 0;
0095
0096 switch (detId.det()) {
0097 case DetId::Tracker:
0098 layer = tTopo->layer(detId);
0099 break;
0100
0101 case DetId::Muon:
0102 switch (detId.subdetId()) {
0103 case MuonSubdetId::DT:
0104 det = TrackQuality::Layer::MuonDT;
0105 layer = DTLayerId(detId).layer();
0106 break;
0107
0108 case MuonSubdetId::CSC:
0109 det = TrackQuality::Layer::MuonCSC;
0110 layer = CSCDetId(detId).layer();
0111 break;
0112
0113 case MuonSubdetId::RPC:
0114 if (RPCDetId(detId).region())
0115 det = TrackQuality::Layer::MuonRPCEndcap;
0116 else
0117 det = TrackQuality::Layer::MuonRPCBarrel;
0118 layer = RPCDetId(detId).layer();
0119 break;
0120
0121 default:
0122
0123 ;
0124 }
0125 break;
0126
0127 default:
0128
0129 ;
0130 }
0131
0132 return DetLayer(det, layer);
0133 }
0134
0135 TrackQuality::TrackQuality(const edm::ParameterSet &config, edm::ConsumesCollector &iC)
0136 : trackerHitAssociatorConfig_(config.getParameter<edm::ParameterSet>("hitAssociator"), std::move(iC)) {}
0137
0138 void TrackQuality::newEvent(const edm::Event &ev, const edm::EventSetup &es) {
0139 associator_ = std::make_unique<TrackerHitAssociator>(ev, trackerHitAssociatorConfig_);
0140 }
0141
0142 void TrackQuality::fillPSetDescription(edm::ParameterSetDescription &desc) {
0143 TrackerHitAssociator::fillPSetDescription(desc);
0144 }
0145
0146 void TrackQuality::evaluate(SimParticleTrail const &spt, reco::TrackBaseRef const &tr, const TrackerTopology *tTopo) {
0147 std::vector<MatchedHit> matchedHits;
0148
0149
0150 for (trackingRecHit_iterator hit = tr->recHitsBegin(); hit != tr->recHitsEnd(); ++hit) {
0151
0152 DetId detId = (*hit)->geographicalId();
0153
0154
0155
0156
0157 if (!(*hit)->isValid()) {
0158 MatchedHit matchedHit;
0159 matchedHit.detId = detId;
0160 matchedHit.simTrackId = NonMatchedTrackId;
0161
0162 switch ((*hit)->getType()) {
0163 case TrackingRecHit::inactive:
0164 matchedHit.state = Layer::Dead;
0165 break;
0166
0167 case TrackingRecHit::bad:
0168 matchedHit.state = Layer::Bad;
0169 break;
0170
0171 default:
0172 matchedHit.state = Layer::Missed;
0173 }
0174 matchedHit.recHitId = hit - tr->recHitsBegin();
0175 matchedHits.push_back(matchedHit);
0176 continue;
0177 }
0178
0179
0180 std::vector<SimHitIdpr> simIds = associator_->associateHitId(**hit);
0181
0182
0183 if (simIds.empty()) {
0184 MatchedHit matchedHit;
0185 matchedHit.detId = detId;
0186 matchedHit.simTrackId = NonMatchedTrackId;
0187 matchedHit.state = Layer::Noise;
0188 matchedHit.recHitId = hit - tr->recHitsBegin();
0189 matchedHits.push_back(matchedHit);
0190 continue;
0191 }
0192
0193
0194 for (std::vector<SimHitIdpr>::const_iterator i = simIds.begin(); i != simIds.end(); ++i) {
0195 MatchedHit matchedHit;
0196 matchedHit.detId = detId;
0197 matchedHit.simTrackId = i->first;
0198 matchedHit.collision = i->second;
0199
0200 if (detId.det() == DetId::Muon)
0201 matchedHit.state = Layer::Unknown;
0202 else
0203
0204 matchedHit.state = Layer::Misassoc;
0205 matchedHit.recHitId = hit - tr->recHitsBegin();
0206 matchedHits.push_back(matchedHit);
0207 }
0208 }
0209
0210
0211 std::stable_sort(matchedHits.begin(), matchedHits.end());
0212
0213
0214 std::stable_sort(matchedHits.begin(), matchedHits.end());
0215
0216
0217 typedef std::multimap<DetLayer, const MatchedHit *> LayerHitMap;
0218 LayerHitMap layerHitMap;
0219
0220
0221 for (std::vector<MatchedHit>::const_iterator hit = matchedHits.begin(); hit != matchedHits.end();) {
0222
0223 const MatchedHit *best = nullptr;
0224
0225
0226 do {
0227
0228 if (!best || statePriorities[hit->state] > statePriorities[best->state] ||
0229 best->simTrackId == NonMatchedTrackId) {
0230 best = &*hit;
0231 }
0232 ++hit;
0233 } while (hit != matchedHits.end() && hit->detId == best->detId);
0234
0235
0236 if (best->simTrackId != NonMatchedTrackId || best->state != Layer::Missed) {
0237 layerHitMap.insert(std::make_pair(getDetLayer(best->detId, tTopo), best));
0238 }
0239 }
0240
0241 layers_.clear();
0242
0243 #ifdef DEBUG_TRACK_QUALITY
0244 std::cout << "---------------------" << std::endl;
0245 #endif
0246
0247 for (LayerHitMap::const_iterator hit = layerHitMap.begin(); hit != layerHitMap.end(); ++hit) {
0248 #ifdef DEBUG_TRACK_QUALITY
0249 std::cout << "detLayer (" << hit->first.first << ", " << hit->first.second << ")"
0250 << " [" << (uint32_t)hit->second->detId << "] sim(" << (int)hit->second->simTrackId << ")"
0251 << " hit(" << hit->second->recHitId << ") -> " << hit->second->state << std::endl;
0252 #endif
0253
0254
0255 Layer *layer = layers_.empty() ? nullptr : &layers_.back();
0256 if (!layer || hit->first.first != layer->subDet || hit->first.second != layer->layer) {
0257 Layer newLayer;
0258 newLayer.subDet = hit->first.first;
0259 newLayer.layer = hit->first.second;
0260 layers_.push_back(newLayer);
0261 layer = &layers_.back();
0262 }
0263
0264
0265 Layer::Hit newHit;
0266 newHit.recHitId = hit->second->recHitId;
0267 newHit.state = hit->second->state;
0268
0269 layer->hits.push_back(newHit);
0270 }
0271 }