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::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   // iterate over reconstructed hits
0150   for (trackingRecHit_iterator hit = tr->recHitsBegin(); hit != tr->recHitsEnd(); ++hit) {
0151     // on which module the hit lies
0152     DetId detId = (*hit)->geographicalId();
0153
0154     // FIXME: check for double-sided modules?
0155
0156     // didn't find a hit on that module
0157     if (!(*hit)->isValid()) {
0158       MatchedHit matchedHit;
0159       matchedHit.detId = detId;
0160       matchedHit.simTrackId = NonMatchedTrackId;
0161       // check why hit wasn't valid and propagate information
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     // find simulated tracks causing hit that was reconstructed
0180     std::vector<SimHitIdpr> simIds = associator_->associateHitId(**hit);
0181
0182     // must be noise or so
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     // register all simulated tracks contributing
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       // RecHit <-> SimHit matcher currently doesn't support muon system
0200       if (detId.det() == DetId::Muon)
0201         matchedHit.state = Layer::Unknown;
0202       else
0203         // assume hit was mismatched (until possible confirmation)
0204         matchedHit.state = Layer::Misassoc;
0205       matchedHit.recHitId = hit - tr->recHitsBegin();
0206       matchedHits.push_back(matchedHit);
0207     }
0208   }
0209
0210   // sort hits found so far by module id
0211   std::stable_sort(matchedHits.begin(), matchedHits.end());
0212
0213   // in case we added missed modules, re-sort
0214   std::stable_sort(matchedHits.begin(), matchedHits.end());
0215
0216   // prepare for ordering results by layer enum and layer/disk number
0217   typedef std::multimap<DetLayer, const MatchedHit *> LayerHitMap;
0218   LayerHitMap layerHitMap;
0219
0220   // iterate over all simulated/reconstructed hits again
0221   for (std::vector<MatchedHit>::const_iterator hit = matchedHits.begin(); hit != matchedHits.end();) {
0222     // we can have multiple reco-to-sim matches per module, find best one
0223     const MatchedHit *best = nullptr;
0224
0225     // this loop iterates over all subsequent hits in the same module
0226     do {
0227       // update our best hit pointer
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     // ignore hit in case track reco was looking at the wrong module
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   // now prepare final collection
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     // find out if we need to start a new layer
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     // create hit and add it to layer
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 }