TrackHistory/src/TrackQuality.cc

0001 /*
0002  *  TrackQuality.cc
0003  *
0004  *  Created by Christophe Saout on 9/25/08.
0005  *  2007 __MyCompanyName__.
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 // #define DEBUG_TRACK_QUALITY
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 static bool operator < (const MatchedHit &hit, DetId detId)
0069 {
0070   return hit.detId < detId;
0071 }
0072 static bool operator < (DetId detId, const MatchedHit &hit)
0073 {
0074   return detId < hit.detId;
0075 }
0076 */
0077 }  // namespace
0078
0079 typedef std::pair<TrackQuality::Layer::SubDet, short int> DetLayer;
0080
0081 // in case multiple hits were found, figure out the highest priority
0082 static const int statePriorities[] = {
0083     /* Unknown */ 3,
0084     /* Good */ 5,
0085     /* Missed */ 0,
0086     /* Noise */ 7,
0087     /* Bad */ 2,
0088     /* Dead */ 4,
0089     /* Shared */ 6,
0090     /* Misassoc */ 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             /* should not get here */
0123             ;
0124       }
0125       break;
0126
0127     default:
0128         /* should not get here */
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::evaluate(SimParticleTrail const &spt, reco::TrackBaseRef const &tr, const TrackerTopology *tTopo) {
0143   std::vector<MatchedHit> matchedHits;
0144
0145   // iterate over reconstructed hits
0146   for (trackingRecHit_iterator hit = tr->recHitsBegin(); hit != tr->recHitsEnd(); ++hit) {
0147     // on which module the hit lies
0148     DetId detId = (*hit)->geographicalId();
0149
0150     // FIXME: check for double-sided modules?
0151
0152     // didn't find a hit on that module
0153     if (!(*hit)->isValid()) {
0154       MatchedHit matchedHit;
0155       matchedHit.detId = detId;
0156       matchedHit.simTrackId = NonMatchedTrackId;
0157       // check why hit wasn't valid and propagate information
0158       switch ((*hit)->getType()) {
0159         case TrackingRecHit::inactive:
0160           matchedHit.state = Layer::Dead;
0161           break;
0162
0163         case TrackingRecHit::bad:
0164           matchedHit.state = Layer::Bad;
0165           break;
0166
0167         default:
0168           matchedHit.state = Layer::Missed;
0169       }
0170       matchedHit.recHitId = hit - tr->recHitsBegin();
0171       matchedHits.push_back(matchedHit);
0172       continue;
0173     }
0174
0175     // find simulated tracks causing hit that was reconstructed
0176     std::vector<SimHitIdpr> simIds = associator_->associateHitId(**hit);
0177
0178     // must be noise or so
0179     if (simIds.empty()) {
0180       MatchedHit matchedHit;
0181       matchedHit.detId = detId;
0182       matchedHit.simTrackId = NonMatchedTrackId;
0183       matchedHit.state = Layer::Noise;
0184       matchedHit.recHitId = hit - tr->recHitsBegin();
0185       matchedHits.push_back(matchedHit);
0186       continue;
0187     }
0188
0189     // register all simulated tracks contributing
0190     for (std::vector<SimHitIdpr>::const_iterator i = simIds.begin(); i != simIds.end(); ++i) {
0191       MatchedHit matchedHit;
0192       matchedHit.detId = detId;
0193       matchedHit.simTrackId = i->first;
0194       matchedHit.collision = i->second;
0195       // RecHit <-> SimHit matcher currently doesn't support muon system
0196       if (detId.det() == DetId::Muon)
0197         matchedHit.state = Layer::Unknown;
0198       else
0199         // assume hit was mismatched (until possible confirmation)
0200         matchedHit.state = Layer::Misassoc;
0201       matchedHit.recHitId = hit - tr->recHitsBegin();
0202       matchedHits.push_back(matchedHit);
0203     }
0204   }
0205
0206   // sort hits found so far by module id
0207   std::stable_sort(matchedHits.begin(), matchedHits.end());
0208
0209   // in case we added missed modules, re-sort
0210   std::stable_sort(matchedHits.begin(), matchedHits.end());
0211
0212   // prepare for ordering results by layer enum and layer/disk number
0213   typedef std::multimap<DetLayer, const MatchedHit *> LayerHitMap;
0214   LayerHitMap layerHitMap;
0215
0216   // iterate over all simulated/reconstructed hits again
0217   for (std::vector<MatchedHit>::const_iterator hit = matchedHits.begin(); hit != matchedHits.end();) {
0218     // we can have multiple reco-to-sim matches per module, find best one
0219     const MatchedHit *best = nullptr;
0220
0221     // this loop iterates over all subsequent hits in the same module
0222     do {
0223       // update our best hit pointer
0224       if (!best || statePriorities[hit->state] > statePriorities[best->state] ||
0225           best->simTrackId == NonMatchedTrackId) {
0226         best = &*hit;
0227       }
0228       ++hit;
0229     } while (hit != matchedHits.end() && hit->detId == best->detId);
0230
0231     // ignore hit in case track reco was looking at the wrong module
0232     if (best->simTrackId != NonMatchedTrackId || best->state != Layer::Missed) {
0233       layerHitMap.insert(std::make_pair(getDetLayer(best->detId, tTopo), best));
0234     }
0235   }
0236
0237   layers_.clear();
0238
0239 #ifdef DEBUG_TRACK_QUALITY
0240   std::cout << "---------------------" << std::endl;
0241 #endif
0242   // now prepare final collection
0243   for (LayerHitMap::const_iterator hit = layerHitMap.begin(); hit != layerHitMap.end(); ++hit) {
0244 #ifdef DEBUG_TRACK_QUALITY
0245     std::cout << "detLayer (" << hit->first.first << ", " << hit->first.second << ")"
0246               << " [" << (uint32_t)hit->second->detId << "] sim(" << (int)hit->second->simTrackId << ")"
0247               << " hit(" << hit->second->recHitId << ") -> " << hit->second->state << std::endl;
0248 #endif
0249
0250     // find out if we need to start a new layer
0251     Layer *layer = layers_.empty() ? nullptr : &layers_.back();
0252     if (!layer || hit->first.first != layer->subDet || hit->first.second != layer->layer) {
0253       Layer newLayer;
0254       newLayer.subDet = hit->first.first;
0255       newLayer.layer = hit->first.second;
0256       layers_.push_back(newLayer);
0257       layer = &layers_.back();
0258     }
0259
0260     // create hit and add it to layer
0261     Layer::Hit newHit;
0262     newHit.recHitId = hit->second->recHitId;
0263     newHit.state = hit->second->state;
0264
0265     layer->hits.push_back(newHit);
0266   }
0267 }