CMS 3D CMS Logo

TrackQuality.cc
Go to the documentation of this file.
1 /*
2  * TrackQuality.cc
3  *
4  * Created by Christophe Saout on 9/25/08.
5  * 2007 __MyCompanyName__.
6  *
7  */
8 
9 #include <algorithm>
10 #include <memory>
11 
12 #include <vector>
13 
15 
18 
28 
33 
35 
37 
38 // #define DEBUG_TRACK_QUALITY
39 
40 namespace {
41  const uint32_t NonMatchedTrackId = (uint32_t)-1;
42 
43  struct MatchedHit {
44  DetId detId;
45  uint32_t simTrackId;
46  EncodedEventId collision;
47  int recHitId;
49 
50  bool operator<(const MatchedHit &other) const {
51  if (detId < other.detId)
52  return true;
53  else if (detId > other.detId)
54  return false;
55  else if (collision < other.collision)
56  return true;
57  else if (other.collision < collision)
58  return false;
59  else
60  return simTrackId < other.simTrackId;
61  }
62 
63  bool operator==(const MatchedHit &other) const {
64  return detId == other.detId && collision == other.collision && simTrackId == other.simTrackId;
65  }
66  };
67  /*
68 static bool operator < (const MatchedHit &hit, DetId detId)
69 {
70  return hit.detId < detId;
71 }
72 static bool operator < (DetId detId, const MatchedHit &hit)
73 {
74  return detId < hit.detId;
75 }
76 */
77 } // namespace
78 
79 typedef std::pair<TrackQuality::Layer::SubDet, short int> DetLayer;
80 
81 // in case multiple hits were found, figure out the highest priority
82 static const int statePriorities[] = {
83  /* Unknown */ 3,
84  /* Good */ 5,
85  /* Missed */ 0,
86  /* Noise */ 7,
87  /* Bad */ 2,
88  /* Dead */ 4,
89  /* Shared */ 6,
90  /* Misassoc */ 1};
91 
94  short int layer = 0;
95 
96  switch (detId.det()) {
97  case DetId::Tracker:
98  layer = tTopo->layer(detId);
99  break;
100 
101  case DetId::Muon:
102  switch (detId.subdetId()) {
103  case MuonSubdetId::DT:
105  layer = DTLayerId(detId).layer();
106  break;
107 
108  case MuonSubdetId::CSC:
110  layer = CSCDetId(detId).layer();
111  break;
112 
113  case MuonSubdetId::RPC:
114  if (RPCDetId(detId).region())
116  else
118  layer = RPCDetId(detId).layer();
119  break;
120 
121  default:
122  /* should not get here */
123  ;
124  }
125  break;
126 
127  default:
128  /* should not get here */
129  ;
130  }
131 
132  return DetLayer(det, layer);
133 }
134 
136  : trackerHitAssociatorConfig_(config.getParameter<edm::ParameterSet>("hitAssociator"), std::move(iC)) {}
137 
139  associator_ = std::make_unique<TrackerHitAssociator>(ev, trackerHitAssociatorConfig_);
140 }
141 
143  std::vector<MatchedHit> matchedHits;
144 
145  // iterate over reconstructed hits
146  for (trackingRecHit_iterator hit = tr->recHitsBegin(); hit != tr->recHitsEnd(); ++hit) {
147  // on which module the hit lies
148  DetId detId = (*hit)->geographicalId();
149 
150  // FIXME: check for double-sided modules?
151 
152  // didn't find a hit on that module
153  if (!(*hit)->isValid()) {
154  MatchedHit matchedHit;
155  matchedHit.detId = detId;
156  matchedHit.simTrackId = NonMatchedTrackId;
157  // check why hit wasn't valid and propagate information
158  switch ((*hit)->getType()) {
160  matchedHit.state = Layer::Dead;
161  break;
162 
163  case TrackingRecHit::bad:
164  matchedHit.state = Layer::Bad;
165  break;
166 
167  default:
168  matchedHit.state = Layer::Missed;
169  }
170  matchedHit.recHitId = hit - tr->recHitsBegin();
171  matchedHits.push_back(matchedHit);
172  continue;
173  }
174 
175  // find simulated tracks causing hit that was reconstructed
176  std::vector<SimHitIdpr> simIds = associator_->associateHitId(**hit);
177 
178  // must be noise or so
179  if (simIds.empty()) {
180  MatchedHit matchedHit;
181  matchedHit.detId = detId;
182  matchedHit.simTrackId = NonMatchedTrackId;
183  matchedHit.state = Layer::Noise;
184  matchedHit.recHitId = hit - tr->recHitsBegin();
185  matchedHits.push_back(matchedHit);
186  continue;
187  }
188 
189  // register all simulated tracks contributing
190  for (std::vector<SimHitIdpr>::const_iterator i = simIds.begin(); i != simIds.end(); ++i) {
191  MatchedHit matchedHit;
192  matchedHit.detId = detId;
193  matchedHit.simTrackId = i->first;
194  matchedHit.collision = i->second;
195  // RecHit <-> SimHit matcher currently doesn't support muon system
196  if (detId.det() == DetId::Muon)
197  matchedHit.state = Layer::Unknown;
198  else
199  // assume hit was mismatched (until possible confirmation)
200  matchedHit.state = Layer::Misassoc;
201  matchedHit.recHitId = hit - tr->recHitsBegin();
202  matchedHits.push_back(matchedHit);
203  }
204  }
205 
206  // sort hits found so far by module id
207  std::stable_sort(matchedHits.begin(), matchedHits.end());
208 
209  // in case we added missed modules, re-sort
210  std::stable_sort(matchedHits.begin(), matchedHits.end());
211 
212  // prepare for ordering results by layer enum and layer/disk number
213  typedef std::multimap<DetLayer, const MatchedHit *> LayerHitMap;
214  LayerHitMap layerHitMap;
215 
216  // iterate over all simulated/reconstructed hits again
217  for (std::vector<MatchedHit>::const_iterator hit = matchedHits.begin(); hit != matchedHits.end();) {
218  // we can have multiple reco-to-sim matches per module, find best one
219  const MatchedHit *best = nullptr;
220 
221  // this loop iterates over all subsequent hits in the same module
222  do {
223  // update our best hit pointer
224  if (!best || statePriorities[hit->state] > statePriorities[best->state] ||
225  best->simTrackId == NonMatchedTrackId) {
226  best = &*hit;
227  }
228  ++hit;
229  } while (hit != matchedHits.end() && hit->detId == best->detId);
230 
231  // ignore hit in case track reco was looking at the wrong module
232  if (best->simTrackId != NonMatchedTrackId || best->state != Layer::Missed) {
233  layerHitMap.insert(std::make_pair(getDetLayer(best->detId, tTopo), best));
234  }
235  }
236 
237  layers_.clear();
238 
239 #ifdef DEBUG_TRACK_QUALITY
240  std::cout << "---------------------" << std::endl;
241 #endif
242  // now prepare final collection
243  for (LayerHitMap::const_iterator hit = layerHitMap.begin(); hit != layerHitMap.end(); ++hit) {
244 #ifdef DEBUG_TRACK_QUALITY
245  std::cout << "detLayer (" << hit->first.first << ", " << hit->first.second << ")"
246  << " [" << (uint32_t)hit->second->detId << "] sim(" << (int)hit->second->simTrackId << ")"
247  << " hit(" << hit->second->recHitId << ") -> " << hit->second->state << std::endl;
248 #endif
249 
250  // find out if we need to start a new layer
251  Layer *layer = layers_.empty() ? nullptr : &layers_.back();
252  if (!layer || hit->first.first != layer->subDet || hit->first.second != layer->layer) {
253  Layer newLayer;
254  newLayer.subDet = hit->first.first;
255  newLayer.layer = hit->first.second;
256  layers_.push_back(newLayer);
257  layer = &layers_.back();
258  }
259 
260  // create hit and add it to layer
261  Layer::Hit newHit;
262  newHit.recHitId = hit->second->recHitId;
263  newHit.state = hit->second->state;
264 
265  layer->hits.push_back(newHit);
266  }
267 }
std::vector< TrackingParticleRef > SimParticleTrail
Definition: TrackQuality.h:30
bool operator<(IOVSyncValue const &iLHS, IOVSyncValue const &iRHS)
Definition: IOVSyncValue.h:34
std::unique_ptr< TrackerHitAssociator > associator_
Definition: TrackQuality.h:81
std::pair< TrackQuality::Layer::SubDet, short int > DetLayer
Definition: TrackQuality.cc:79
TrackerHitAssociator::Config trackerHitAssociatorConfig_
Definition: TrackQuality.h:80
Definition: config.py:1
int layer() const
Definition: CSCDetId.h:56
unsigned int layer(const DetId &id) const
void evaluate(SimParticleTrail const &, reco::TrackBaseRef const &, const TrackerTopology *tTopo)
Compute information about the track reconstruction quality.
bool operator==(const QGLikelihoodParameters &lhs, const QGLikelihoodCategory &rhs)
Test if parameters are compatible with category.
trackingRecHit_iterator recHitsEnd() const
Iterator to last hit on the track.
Definition: Track.h:91
DetLayer getDetLayer(DetId detId, const TrackerTopology *tTopo)
Definition: TrackQuality.cc:92
std::vector< Layer > layers_
Definition: TrackQuality.h:83
trackingRecHit_iterator recHitsBegin() const
Iterator to first hit on the track.
Definition: Track.h:88
Definition: DetId.h:17
static constexpr int RPC
Definition: MuonSubdetId.h:13
const Layer & layer(unsigned int index) const
Return information about the given layer by index.
Definition: TrackQuality.h:77
int layer() const
Return the layer number.
Definition: DTLayerId.h:45
HLT enums.
std::vector< Hit > hits
Definition: TrackQuality.h:56
void newEvent(const edm::Event &, const edm::EventSetup &)
Pre-process event information (for accessing reconstruction information)
int layer() const
Definition: RPCDetId.h:85
static constexpr int DT
Definition: MuonSubdetId.h:11
TrackQuality(const edm::ParameterSet &, edm::ConsumesCollector &iC)
Constructor by pset.
static constexpr int CSC
Definition: MuonSubdetId.h:12
def move(src, dest)
Definition: eostools.py:511
static const int statePriorities[]
Definition: TrackQuality.cc:82