CMS 3D CMS Logo

MuonAssociatorByHitsHelper.cc
Go to the documentation of this file.
12 #include <sstream>
13 
14 using namespace reco;
15 using namespace std;
16 
18  : includeZeroHitMuons(conf.getParameter<bool>("includeZeroHitMuons")),
19  acceptOneStubMatchings(conf.getParameter<bool>("acceptOneStubMatchings")),
20  UseTracker(conf.getParameter<bool>("UseTracker")),
21  UseMuon(conf.getParameter<bool>("UseMuon")),
22  AbsoluteNumberOfHits_track(conf.getParameter<bool>("AbsoluteNumberOfHits_track")),
23  NHitCut_track(conf.getParameter<unsigned int>("NHitCut_track")),
24  EfficiencyCut_track(conf.getParameter<double>("EfficiencyCut_track")),
25  PurityCut_track(conf.getParameter<double>("PurityCut_track")),
26  AbsoluteNumberOfHits_muon(conf.getParameter<bool>("AbsoluteNumberOfHits_muon")),
27  NHitCut_muon(conf.getParameter<unsigned int>("NHitCut_muon")),
28  EfficiencyCut_muon(conf.getParameter<double>("EfficiencyCut_muon")),
29  PurityCut_muon(conf.getParameter<double>("PurityCut_muon")),
30  UsePixels(conf.getParameter<bool>("UsePixels")),
31  UseGrouped(conf.getParameter<bool>("UseGrouped")),
32  UseSplitting(conf.getParameter<bool>("UseSplitting")),
33  ThreeHitTracksAreSpecial(conf.getParameter<bool>("ThreeHitTracksAreSpecial")),
34  dumpDT(conf.getParameter<bool>("dumpDT")) {
35  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "constructing MuonAssociatorByHitsHelper" << conf.dump();
36 
37  // up to the user in the other cases - print a message
38  if (UseTracker)
39  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "\n UseTracker = TRUE : Tracker SimHits and RecHits WILL be "
40  "counted";
41  else
42  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "\n UseTracker = FALSE : Tracker SimHits and RecHits WILL NOT be "
43  "counted";
44 
45  // up to the user in the other cases - print a message
46  if (UseMuon)
47  edm::LogVerbatim("MuonAssociatorByHitsHelper") << " UseMuon = TRUE : Muon SimHits and RecHits WILL be counted";
48  else
49  edm::LogVerbatim("MuonAssociatorByHitsHelper")
50  << " UseMuon = FALSE : Muon SimHits and RecHits WILL NOT be counted" << endl;
51 
52  // check consistency of the configuration when allowing zero-hit muon matching
53  // (counting invalid hits)
54  if (includeZeroHitMuons) {
55  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "\n includeZeroHitMuons = TRUE"
56  << "\n ==> (re)set NHitCut_muon = 0, PurityCut_muon = 0, "
57  "EfficiencyCut_muon = 0"
58  << endl;
59  NHitCut_muon = 0;
60  PurityCut_muon = 0.;
61  EfficiencyCut_muon = 0.;
62  }
63 }
64 
66  const TrackHitsCollection &tC,
67  const edm::RefVector<TrackingParticleCollection> &TPCollectionH,
68  const Resources &resources) const {
69  auto tTopo = resources.tTopo_;
70  auto trackertruth = resources.trackerHitAssoc_;
71  auto const &csctruth = *resources.cscHitAssoc_;
72  auto const &dttruth = *resources.dtHitAssoc_;
73  auto const &rpctruth = *resources.rpcHitAssoc_;
74  auto const &gemtruth = *resources.gemHitAssoc_;
75 
76  int tracker_nshared = 0;
77  int muon_nshared = 0;
78  int global_nshared = 0;
79 
80  double tracker_quality = 0;
81  double tracker_quality_cut;
83  tracker_quality_cut = static_cast<double>(NHitCut_track);
84  else
85  tracker_quality_cut = PurityCut_track;
86 
87  double muon_quality = 0;
88  double muon_quality_cut;
90  muon_quality_cut = static_cast<double>(NHitCut_muon);
91  else
92  muon_quality_cut = PurityCut_muon;
93 
94  double global_quality = 0;
95 
96  MapOfMatchedIds tracker_matchedIds_valid, muon_matchedIds_valid;
97  MapOfMatchedIds tracker_matchedIds_INVALID, muon_matchedIds_INVALID;
98 
99  IndexAssociation outputCollection;
100  bool printRtS(true);
101 
103  tPC.reserve(TPCollectionH.size());
104  for (auto const &ref : TPCollectionH) {
105  tPC.push_back(*ref);
106  }
107 
108  if (resources.diagnostics_) {
109  resources.diagnostics_(tC, tPC);
110  }
111 
112  int tindex = 0;
113  for (TrackHitsCollection::const_iterator track = tC.begin(); track != tC.end(); track++, tindex++) {
114  edm::LogVerbatim("MuonAssociatorByHitsHelper")
115  << "\n"
116  << "reco::Track " << tindex << ", number of RecHits = " << (track->second - track->first) << "\n";
117  tracker_matchedIds_valid.clear();
118  muon_matchedIds_valid.clear();
119 
120  tracker_matchedIds_INVALID.clear();
121  muon_matchedIds_INVALID.clear();
122 
123  bool this_track_matched = false;
124  int n_matching_simhits = 0;
125 
126  // all hits = valid +INVALID
127  int n_all = 0;
128  int n_tracker_all = 0;
129  int n_dt_all = 0;
130  int n_csc_all = 0;
131  int n_rpc_all = 0;
132  int n_gem_all = 0;
133 
134  int n_valid = 0;
135  int n_tracker_valid = 0;
136  int n_muon_valid = 0;
137  int n_dt_valid = 0;
138  int n_csc_valid = 0;
139  int n_rpc_valid = 0;
140  int n_gem_valid = 0;
141 
142  int n_tracker_matched_valid = 0;
143  int n_muon_matched_valid = 0;
144  int n_dt_matched_valid = 0;
145  int n_csc_matched_valid = 0;
146  int n_rpc_matched_valid = 0;
147  int n_gem_matched_valid = 0;
148 
149  int n_INVALID = 0;
150  int n_tracker_INVALID = 0;
151  int n_muon_INVALID = 0;
152  int n_dt_INVALID = 0;
153  int n_csc_INVALID = 0;
154  int n_rpc_INVALID = 0;
155  int n_gem_INVALID = 0;
156 
157  int n_tracker_matched_INVALID = 0;
158  int n_muon_matched_INVALID = 0;
159  int n_dt_matched_INVALID = 0;
160  int n_csc_matched_INVALID = 0;
161  int n_rpc_matched_INVALID = 0;
162  int n_gem_matched_INVALID = 0;
163 
164  printRtS = true;
165  getMatchedIds(tracker_matchedIds_valid,
166  muon_matchedIds_valid,
167  tracker_matchedIds_INVALID,
168  muon_matchedIds_INVALID,
169  n_tracker_valid,
170  n_dt_valid,
171  n_csc_valid,
172  n_rpc_valid,
173  n_gem_valid,
174  n_tracker_matched_valid,
175  n_dt_matched_valid,
176  n_csc_matched_valid,
177  n_rpc_matched_valid,
178  n_gem_matched_valid,
179  n_tracker_INVALID,
180  n_dt_INVALID,
181  n_csc_INVALID,
182  n_rpc_INVALID,
183  n_gem_INVALID,
184  n_tracker_matched_INVALID,
185  n_dt_matched_INVALID,
186  n_csc_matched_INVALID,
187  n_rpc_matched_INVALID,
188  n_gem_matched_INVALID,
189  track->first,
190  track->second,
191  trackertruth,
192  dttruth,
193  csctruth,
194  rpctruth,
195  gemtruth,
196  printRtS,
197  tTopo);
198 
199  n_matching_simhits = tracker_matchedIds_valid.size() + muon_matchedIds_valid.size() +
200  tracker_matchedIds_INVALID.size() + muon_matchedIds_INVALID.size();
201 
202  n_muon_valid = n_dt_valid + n_csc_valid + n_rpc_valid + n_gem_valid;
203  n_valid = n_tracker_valid + n_muon_valid;
204  n_muon_INVALID = n_dt_INVALID + n_csc_INVALID + n_rpc_INVALID + n_gem_INVALID;
205  n_INVALID = n_tracker_INVALID + n_muon_INVALID;
206 
207  // all used hits (valid+INVALID), defined by UseTracker, UseMuon
208  n_tracker_all = n_tracker_valid + n_tracker_INVALID;
209  n_dt_all = n_dt_valid + n_dt_INVALID;
210  n_csc_all = n_csc_valid + n_csc_INVALID;
211  n_rpc_all = n_rpc_valid + n_rpc_INVALID;
212  n_gem_all = n_gem_valid + n_gem_INVALID;
213  n_all = n_valid + n_INVALID;
214 
215  n_muon_matched_valid = n_dt_matched_valid + n_csc_matched_valid + n_rpc_matched_valid + n_gem_matched_valid;
216  n_muon_matched_INVALID =
217  n_dt_matched_INVALID + n_csc_matched_INVALID + n_rpc_matched_INVALID + n_gem_matched_INVALID;
218 
219  // selected hits are set initially to valid hits
220  int n_tracker_selected_hits = n_tracker_valid;
221  int n_muon_selected_hits = n_muon_valid;
222  int n_dt_selected_hits = n_dt_valid;
223  int n_csc_selected_hits = n_csc_valid;
224  int n_rpc_selected_hits = n_rpc_valid;
225  int n_gem_selected_hits = n_gem_valid;
226 
227  // matched hits are a subsample of the selected hits
228  int n_tracker_matched = n_tracker_matched_valid;
229  int n_muon_matched = n_muon_matched_valid;
230  int n_dt_matched = n_dt_matched_valid;
231  int n_csc_matched = n_csc_matched_valid;
232  int n_rpc_matched = n_rpc_matched_valid;
233  int n_gem_matched = n_gem_matched_valid;
234 
235  std::string InvMuonHits, ZeroHitMuon;
236 
237  if (includeZeroHitMuons && n_muon_valid == 0 && n_muon_INVALID != 0) {
238  // selected muon hits = INVALID when (useZeroHitMuons == True) and track
239  // has no valid muon hits
240 
241  InvMuonHits = " ***INVALID MUON HITS***";
242  ZeroHitMuon = " ***ZERO-HIT MUON***";
243 
244  n_muon_selected_hits = n_muon_INVALID;
245  n_dt_selected_hits = n_dt_INVALID;
246  n_csc_selected_hits = n_csc_INVALID;
247  n_rpc_selected_hits = n_rpc_INVALID;
248  n_gem_selected_hits = n_gem_INVALID;
249 
250  n_muon_matched = n_muon_matched_INVALID;
251  n_dt_matched = n_dt_matched_INVALID;
252  n_csc_matched = n_csc_matched_INVALID;
253  n_rpc_matched = n_rpc_matched_INVALID;
254  n_gem_matched = n_gem_matched_INVALID;
255  }
256 
257  int n_selected_hits = n_tracker_selected_hits + n_muon_selected_hits;
258  int n_matched = n_tracker_matched + n_muon_matched;
259 
260  edm::LogVerbatim("MuonAssociatorByHitsHelper")
261  << "\n"
262  << "# TrackingRecHits: " << (track->second - track->first) << "\n"
263  << "# used RecHits = " << n_all << " (" << n_tracker_all << "/" << n_dt_all << "/" << n_csc_all << "/"
264  << n_rpc_all << "/" << n_gem_all << " in Tracker/DT/CSC/RPC/GEM)"
265  << ", obtained from " << n_matching_simhits << " SimHits"
266  << "\n"
267  << "# selected RecHits = " << n_selected_hits << " (" << n_tracker_selected_hits << "/" << n_dt_selected_hits
268  << "/" << n_csc_selected_hits << "/" << n_rpc_selected_hits << "/" << n_gem_selected_hits
269  << " in Tracker/DT/CSC/RPC/GEM)" << InvMuonHits << "\n"
270  << "# matched RecHits = " << n_matched << " (" << n_tracker_matched << "/" << n_dt_matched << "/"
271  << n_csc_matched << "/" << n_rpc_matched << "/" << n_gem_matched << " in Tracker/DT/CSC/RPC/GEM)";
272 
273  if (n_all > 0 && n_matching_simhits == 0)
274  edm::LogVerbatim("MuonAssociatorByHitsHelper")
275  << "*** WARNING in MuonAssociatorByHitsHelper::associateRecoToSim: "
276  "no matching PSimHit found for this reco::Track !";
277 
278  if (n_matching_simhits != 0) {
279  edm::LogVerbatim("MuonAssociatorByHitsHelper")
280  << "\n"
281  << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
282  "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
283  << "\n"
284  << "reco::Track " << tindex << ZeroHitMuon << "\n\t"
285  << "made of " << n_selected_hits << " selected RecHits (tracker:" << n_tracker_selected_hits
286  << "/muons:" << n_muon_selected_hits << ")";
287 
288  int tpindex = 0;
289  for (TrackingParticleCollection::const_iterator trpart = tPC.begin(); trpart != tPC.end(); ++trpart, ++tpindex) {
290  tracker_nshared = getShared(tracker_matchedIds_valid, trpart);
291  muon_nshared = getShared(muon_matchedIds_valid, trpart);
292 
293  if (includeZeroHitMuons && n_muon_valid == 0 && n_muon_INVALID != 0)
294  muon_nshared = getShared(muon_matchedIds_INVALID, trpart);
295 
296  global_nshared = tracker_nshared + muon_nshared;
297 
299  tracker_quality = static_cast<double>(tracker_nshared);
300  else if (n_tracker_selected_hits != 0)
301  tracker_quality = (static_cast<double>(tracker_nshared) / static_cast<double>(n_tracker_selected_hits));
302  else
303  tracker_quality = 0;
304 
306  muon_quality = static_cast<double>(muon_nshared);
307  else if (n_muon_selected_hits != 0)
308  muon_quality = (static_cast<double>(muon_nshared) / static_cast<double>(n_muon_selected_hits));
309  else
310  muon_quality = 0;
311 
312  // global_quality used to order the matching TPs
313  if (n_selected_hits != 0) {
315  global_quality = global_nshared;
316  else
317  global_quality = (static_cast<double>(global_nshared) / static_cast<double>(n_selected_hits));
318  } else
319  global_quality = 0;
320 
321  bool trackerOk = false;
322  if (n_tracker_selected_hits != 0) {
323  if (tracker_quality > tracker_quality_cut)
324  trackerOk = true;
325  // if a track has just 3 hits in the Tracker we require that all 3
326  // hits are shared
327  if (ThreeHitTracksAreSpecial && n_tracker_selected_hits == 3 && tracker_nshared < 3)
328  trackerOk = false;
329  }
330 
331  bool muonOk = false;
332  if (n_muon_selected_hits != 0) {
333  if (muon_quality > muon_quality_cut)
334  muonOk = true;
335  }
336 
337  // (matchOk) has to account for different track types (tracker-only,
338  // standalone muons, global muons)
339  bool matchOk = trackerOk || muonOk;
340 
341  // only for global muons: match both tracker and muon stub unless
342  // (acceptOneStubMatchings==true)
343  if (!acceptOneStubMatchings && n_tracker_selected_hits != 0 && n_muon_selected_hits != 0)
344  matchOk = trackerOk && muonOk;
345 
346  if (matchOk) {
347  outputCollection[tindex].push_back(IndexMatch(tpindex, global_quality));
348  this_track_matched = true;
349 
350  edm::LogVerbatim("MuonAssociatorByHitsHelper")
351  << "\n\t"
352  << " **MATCHED** with quality = " << global_quality << " (tracker: " << tracker_quality
353  << " / muon: " << muon_quality << ")"
354  << "\n\t"
355  << " N shared hits = " << global_nshared << " (tracker: " << tracker_nshared
356  << " / muon: " << muon_nshared << ")"
357  << "\n"
358  << " to: TrackingParticle " << tpindex << ", q = " << (*trpart).charge() << ", p = " << (*trpart).p()
359  << ", pT = " << (*trpart).pt() << ", eta = " << (*trpart).eta() << ", phi = " << (*trpart).phi() << "\n\t"
360  << " pdg code = " << (*trpart).pdgId() << ", made of " << (*trpart).numberOfHits() << " PSimHits"
361  << " from " << (*trpart).g4Tracks().size() << " SimTrack:";
362  for (TrackingParticle::g4t_iterator g4T = (*trpart).g4Track_begin(); g4T != (*trpart).g4Track_end(); ++g4T) {
363  edm::LogVerbatim("MuonAssociatorByHitsHelper")
364  << "\t"
365  << " Id:" << (*g4T).trackId() << "/Evt:(" << (*g4T).eventId().event() << ","
366  << (*g4T).eventId().bunchCrossing() << ")";
367  }
368  } else {
369  // print something only if this TrackingParticle shares some hits with
370  // the current reco::Track
371  if (global_nshared != 0)
372  edm::LogVerbatim("MuonAssociatorByHitsHelper")
373  << "\n\t"
374  << " NOT matched to TrackingParticle " << tpindex << " with quality = " << global_quality
375  << " (tracker: " << tracker_quality << " / muon: " << muon_quality << ")"
376  << "\n"
377  << " N shared hits = " << global_nshared << " (tracker: " << tracker_nshared
378  << " / muon: " << muon_nshared << ")";
379  }
380 
381  } // loop over TrackingParticle
382 
383  if (!this_track_matched) {
384  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "\n"
385  << " NOT matched to any TrackingParticle";
386  }
387 
388  edm::LogVerbatim("MuonAssociatorByHitsHelper")
389  << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
390  "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
391  << "\n";
392 
393  } // if(n_matching_simhits != 0)
394 
395  } // loop over reco::Track
396 
397  if (tC.empty())
398  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "0 reconstructed tracks (-->> 0 associated !)";
399 
400  for (IndexAssociation::iterator it = outputCollection.begin(), ed = outputCollection.end(); it != ed; ++it) {
401  std::sort(it->second.begin(), it->second.end());
402  }
403  return outputCollection;
404 }
405 
407  const TrackHitsCollection &tC,
408  const edm::RefVector<TrackingParticleCollection> &TPCollectionH,
409  const Resources &resources) const {
410  auto tTopo = resources.tTopo_;
411  auto trackertruth = resources.trackerHitAssoc_;
412  auto const &csctruth = *resources.cscHitAssoc_;
413  auto const &dttruth = *resources.dtHitAssoc_;
414  auto const &rpctruth = *resources.rpcHitAssoc_;
415  auto const &gemtruth = *resources.gemHitAssoc_;
416 
417  int tracker_nshared = 0;
418  int muon_nshared = 0;
419  int global_nshared = 0;
420 
421  double tracker_quality = 0;
422  double tracker_quality_cut;
424  tracker_quality_cut = static_cast<double>(NHitCut_track);
425  else
426  tracker_quality_cut = EfficiencyCut_track;
427 
428  double muon_quality = 0;
429  double muon_quality_cut;
431  muon_quality_cut = static_cast<double>(NHitCut_muon);
432  else
433  muon_quality_cut = EfficiencyCut_muon;
434 
435  double global_quality = 0;
436 
437  double tracker_purity = 0;
438  double muon_purity = 0;
439  double global_purity = 0;
440 
441  MapOfMatchedIds tracker_matchedIds_valid, muon_matchedIds_valid;
442  MapOfMatchedIds tracker_matchedIds_INVALID, muon_matchedIds_INVALID;
443 
444  IndexAssociation outputCollection;
445 
446  bool printRtS(true);
447 
449  tPC.reserve(TPCollectionH.size());
450  for (auto const &ref : TPCollectionH) {
451  tPC.push_back(*ref);
452  }
453 
454  bool any_trackingParticle_matched = false;
455 
456  int tindex = 0;
457  for (TrackHitsCollection::const_iterator track = tC.begin(); track != tC.end(); track++, tindex++) {
458  if (printRtS)
459  edm::LogVerbatim("MuonAssociatorByHitsHelper")
460  << "\n"
461  << "reco::Track " << tindex << ", number of RecHits = " << (track->second - track->first) << "\n";
462 
463  tracker_matchedIds_valid.clear();
464  muon_matchedIds_valid.clear();
465 
466  tracker_matchedIds_INVALID.clear();
467  muon_matchedIds_INVALID.clear();
468 
469  int n_matching_simhits = 0;
470 
471  // all hits = valid +INVALID
472  int n_all = 0;
473  int n_tracker_all = 0;
474  int n_dt_all = 0;
475  int n_csc_all = 0;
476  int n_rpc_all = 0;
477  int n_gem_all = 0;
478 
479  int n_valid = 0;
480  int n_tracker_valid = 0;
481  int n_muon_valid = 0;
482  int n_dt_valid = 0;
483  int n_csc_valid = 0;
484  int n_rpc_valid = 0;
485  int n_gem_valid = 0;
486 
487  int n_tracker_matched_valid = 0;
488  int n_muon_matched_valid = 0;
489  int n_dt_matched_valid = 0;
490  int n_csc_matched_valid = 0;
491  int n_rpc_matched_valid = 0;
492  int n_gem_matched_valid = 0;
493 
494  int n_INVALID = 0;
495  int n_tracker_INVALID = 0;
496  int n_muon_INVALID = 0;
497  int n_dt_INVALID = 0;
498  int n_csc_INVALID = 0;
499  int n_rpc_INVALID = 0;
500  int n_gem_INVALID = 0;
501 
502  int n_tracker_matched_INVALID = 0;
503  int n_muon_matched_INVALID = 0;
504  int n_dt_matched_INVALID = 0;
505  int n_csc_matched_INVALID = 0;
506  int n_rpc_matched_INVALID = 0;
507  int n_gem_matched_INVALID = 0;
508 
509  printRtS = false;
510  getMatchedIds(tracker_matchedIds_valid,
511  muon_matchedIds_valid,
512  tracker_matchedIds_INVALID,
513  muon_matchedIds_INVALID,
514  n_tracker_valid,
515  n_dt_valid,
516  n_csc_valid,
517  n_rpc_valid,
518  n_gem_valid,
519  n_tracker_matched_valid,
520  n_dt_matched_valid,
521  n_csc_matched_valid,
522  n_rpc_matched_valid,
523  n_gem_matched_valid,
524  n_tracker_INVALID,
525  n_dt_INVALID,
526  n_csc_INVALID,
527  n_rpc_INVALID,
528  n_gem_INVALID,
529  n_tracker_matched_INVALID,
530  n_dt_matched_INVALID,
531  n_csc_matched_INVALID,
532  n_rpc_matched_INVALID,
533  n_gem_matched_INVALID,
534  track->first,
535  track->second,
536  trackertruth,
537  dttruth,
538  csctruth,
539  rpctruth,
540  gemtruth,
541  printRtS,
542  tTopo);
543 
544  n_matching_simhits = tracker_matchedIds_valid.size() + muon_matchedIds_valid.size() +
545  tracker_matchedIds_INVALID.size() + muon_matchedIds_INVALID.size();
546 
547  n_muon_valid = n_dt_valid + n_csc_valid + n_rpc_valid + n_gem_valid;
548  n_valid = n_tracker_valid + n_muon_valid;
549  n_muon_INVALID = n_dt_INVALID + n_csc_INVALID + n_rpc_INVALID + n_gem_INVALID;
550  n_INVALID = n_tracker_INVALID + n_muon_INVALID;
551 
552  // all used hits (valid+INVALID), defined by UseTracker, UseMuon
553  n_tracker_all = n_tracker_valid + n_tracker_INVALID;
554  n_dt_all = n_dt_valid + n_dt_INVALID;
555  n_csc_all = n_csc_valid + n_csc_INVALID;
556  n_rpc_all = n_rpc_valid + n_rpc_INVALID;
557  n_gem_all = n_gem_valid + n_gem_INVALID;
558  n_all = n_valid + n_INVALID;
559 
560  n_muon_matched_valid = n_dt_matched_valid + n_csc_matched_valid + n_rpc_matched_valid + n_gem_matched_valid;
561  n_muon_matched_INVALID =
562  n_dt_matched_INVALID + n_csc_matched_INVALID + n_rpc_matched_INVALID + n_gem_matched_INVALID;
563 
564  // selected hits are set initially to valid hits
565  int n_tracker_selected_hits = n_tracker_valid;
566  int n_muon_selected_hits = n_muon_valid;
567  int n_dt_selected_hits = n_dt_valid;
568  int n_csc_selected_hits = n_csc_valid;
569  int n_rpc_selected_hits = n_rpc_valid;
570  int n_gem_selected_hits = n_gem_valid;
571 
572  // matched hits are a subsample of the selected hits
573  int n_tracker_matched = n_tracker_matched_valid;
574  int n_muon_matched = n_muon_matched_valid;
575  int n_dt_matched = n_dt_matched_valid;
576  int n_csc_matched = n_csc_matched_valid;
577  int n_rpc_matched = n_rpc_matched_valid;
578  int n_gem_matched = n_gem_matched_valid;
579 
580  std::string InvMuonHits, ZeroHitMuon;
581 
582  if (includeZeroHitMuons && n_muon_valid == 0 && n_muon_INVALID != 0) {
583  // selected muon hits = INVALID when (useZeroHitMuons == True) and track
584  // has no valid muon hits
585 
586  InvMuonHits = " ***INVALID MUON HITS***";
587  ZeroHitMuon = " ***ZERO-HIT MUON***";
588 
589  n_muon_selected_hits = n_muon_INVALID;
590  n_dt_selected_hits = n_dt_INVALID;
591  n_csc_selected_hits = n_csc_INVALID;
592  n_rpc_selected_hits = n_rpc_INVALID;
593  n_gem_selected_hits = n_gem_INVALID;
594 
595  n_muon_matched = n_muon_matched_INVALID;
596  n_dt_matched = n_dt_matched_INVALID;
597  n_csc_matched = n_csc_matched_INVALID;
598  n_rpc_matched = n_rpc_matched_INVALID;
599  n_gem_matched = n_gem_matched_INVALID;
600  }
601 
602  int n_selected_hits = n_tracker_selected_hits + n_muon_selected_hits;
603  int n_matched = n_tracker_matched + n_muon_matched;
604 
605  if (printRtS)
606  edm::LogVerbatim("MuonAssociatorByHitsHelper")
607  << "\n"
608  << "# TrackingRecHits: " << (track->second - track->first) << "\n"
609  << "# used RecHits = " << n_all << " (" << n_tracker_all << "/" << n_dt_all << "/" << n_csc_all << "/"
610  << n_rpc_all << "/" << n_gem_all << " in Tracker/DT/CSC/RPC/GEM)"
611  << ", obtained from " << n_matching_simhits << " SimHits"
612  << "\n"
613  << "# selected RecHits = " << n_selected_hits << " (" << n_tracker_selected_hits << "/" << n_dt_selected_hits
614  << "/" << n_csc_selected_hits << "/" << n_rpc_selected_hits << "/" << n_gem_selected_hits
615  << " in Tracker/DT/CSC/RPC/GEM)" << InvMuonHits << "\n"
616  << "# matched RecHits = " << n_matched << " (" << n_tracker_matched << "/" << n_dt_matched << "/"
617  << n_csc_matched << "/" << n_rpc_matched << "/" << n_gem_matched << " in Tracker/DT/CSC/RPC/GEM)";
618 
619  if (printRtS && n_all > 0 && n_matching_simhits == 0)
620  edm::LogVerbatim("MuonAssociatorByHitsHelper")
621  << "*** WARNING in MuonAssociatorByHitsHelper::associateSimToReco: "
622  "no matching PSimHit found for this reco::Track !";
623 
624  if (n_matching_simhits != 0) {
625  int tpindex = 0;
626  for (TrackingParticleCollection::const_iterator trpart = tPC.begin(); trpart != tPC.end(); ++trpart, ++tpindex) {
627  // int n_tracker_simhits = 0;
628  int n_tracker_recounted_simhits = 0;
629  int n_muon_simhits = 0;
630  int n_global_simhits = 0;
631  // std::vector<PSimHit> tphits;
632 
633  int n_tracker_selected_simhits = 0;
634  int n_muon_selected_simhits = 0;
635  int n_global_selected_simhits = 0;
636 
637  // shared hits are counted over the selected ones
638  tracker_nshared = getShared(tracker_matchedIds_valid, trpart);
639  muon_nshared = getShared(muon_matchedIds_valid, trpart);
640 
641  if (includeZeroHitMuons && n_muon_valid == 0 && n_muon_INVALID != 0)
642  muon_nshared = getShared(muon_matchedIds_INVALID, trpart);
643 
644  global_nshared = tracker_nshared + muon_nshared;
645  if (global_nshared == 0)
646  continue; // if this TP shares no hits with the current reco::Track
647  // loop over
648 
649  // This does not work with the new TP interface
650  /*
651  for(std::vector<PSimHit>::const_iterator TPhit =
652  trpart->pSimHit_begin(); TPhit != trpart->pSimHit_end(); TPhit++) {
653  DetId dId = DetId(TPhit->detUnitId());
654  DetId::Detector detector = dId.det();
655 
656  if (detector == DetId::Tracker) {
657  n_tracker_simhits++;
658 
659  unsigned int subdetId = static_cast<unsigned int>(dId.subdetId());
660  if (!UsePixels && (subdetId==PixelSubdetector::PixelBarrel ||
661  subdetId==PixelSubdetector::PixelEndcap) ) continue;
662 
663  SiStripDetId* stripDetId = 0;
664  if (subdetId==SiStripDetId::TIB||subdetId==SiStripDetId::TOB||
665  subdetId==SiStripDetId::TID||subdetId==SiStripDetId::TEC)
666  stripDetId= new SiStripDetId(dId);
667 
668  bool newhit = true;
669  for(std::vector<PSimHit>::const_iterator TPhitOK = tphits.begin();
670  TPhitOK != tphits.end(); TPhitOK++) { DetId dIdOK =
671  DetId(TPhitOK->detUnitId());
672  //no grouped, no splitting
673  if (!UseGrouped && !UseSplitting)
674  if (tTopo->layer(dId)==tTopo->layer(dIdOK) &&
675  dId.subdetId()==dIdOK.subdetId()) newhit = false;
676  //no grouped, splitting
677  if (!UseGrouped && UseSplitting)
678  if (tTopo->layer(dId)==tTopo->layer(dIdOK) &&
679  dId.subdetId()==dIdOK.subdetId() &&
680  (stripDetId==0 ||
681  stripDetId->partnerDetId()!=dIdOK.rawId())) newhit = false;
682  //grouped, no splitting
683  if (UseGrouped && !UseSplitting)
684  if ( tTopo->layer(dId)== tTopo->layer(dIdOK) &&
685  dId.subdetId()==dIdOK.subdetId() &&
686  stripDetId!=0 && stripDetId->partnerDetId()==dIdOK.rawId())
687  newhit = false;
688  //grouped, splitting
689  if (UseGrouped && UseSplitting)
690  newhit = true;
691  }
692  if (newhit) {
693  tphits.push_back(*TPhit);
694  }
695  delete stripDetId;
696  }
697  else if (detector == DetId::Muon) {
698  n_muon_simhits++;
699 
700  // discard BAD CSC chambers (ME4/2) from hit counting
701  if (dId.subdetId() == MuonSubdetId::CSC) {
702  if (csctruth.cscBadChambers->isInBadChamber(CSCDetId(dId))) {
703  // edm::LogVerbatim("MuonAssociatorByHitsHelper")<<"This PSimHit
704  is in a BAD CSC chamber, CSCDetId = "<<CSCDetId(dId); n_muon_simhits--;
705  }
706  }
707 
708  }
709  }
710  */
711  // n_tracker_recounted_simhits = tphits.size();
712 
713  // adapt to new TP interface: this gives the total number of hits in
714  // tracker
715  // should reproduce the behaviour of UseGrouped=UseSplitting=.true.
716  n_tracker_recounted_simhits = trpart->numberOfTrackerHits();
717  // numberOfHits() gives the total number of hits (tracker + muons)
718  n_muon_simhits = trpart->numberOfHits() - trpart->numberOfTrackerHits();
719 
720  // Handle the case of TrackingParticles that don't have PSimHits inside,
721  // e.g. because they were made on RECOSIM only.
722  if (trpart->numberOfHits() == 0) {
723  // FIXME this can be made better, counting the digiSimLinks associated
724  // to this TP, but perhaps it's not worth it
725  // n_tracker_recounted_simhits = tracker_nshared;
726  // n_muon_simhits = muon_nshared;
727  // ---> on RECOSIM when the AbsoluteNumberOfHits_muon=True this always
728  // obtains quality=1,
729  // hence no sorting is possible in case of duplicate matchings
730  // ---> reset these variables to 1 so to keep the number of shared
731  // hits as ranking criterion
732  n_tracker_recounted_simhits = 1;
733  n_muon_simhits = 1;
734  }
735  n_global_simhits = n_tracker_recounted_simhits + n_muon_simhits;
736 
737  if (UseMuon) {
738  n_muon_selected_simhits = n_muon_simhits;
739  n_global_selected_simhits = n_muon_selected_simhits;
740  }
741  if (UseTracker) {
742  n_tracker_selected_simhits = n_tracker_recounted_simhits;
743  n_global_selected_simhits += n_tracker_selected_simhits;
744  }
745 
747  tracker_quality = static_cast<double>(tracker_nshared);
748  else if (n_tracker_selected_simhits != 0)
749  tracker_quality = static_cast<double>(tracker_nshared) / static_cast<double>(n_tracker_selected_simhits);
750  else
751  tracker_quality = 0;
752 
754  muon_quality = static_cast<double>(muon_nshared);
755  else if (n_muon_selected_simhits != 0)
756  muon_quality = static_cast<double>(muon_nshared) / static_cast<double>(n_muon_selected_simhits);
757  else
758  muon_quality = 0;
759 
760  // global_quality used to order the matching tracks
761  if (n_global_selected_simhits != 0) {
763  global_quality = global_nshared;
764  else
765  global_quality = static_cast<double>(global_nshared) / static_cast<double>(n_global_selected_simhits);
766  } else
767  global_quality = 0;
768 
769  // global purity
770  if (n_selected_hits != 0) {
772  global_purity = global_nshared;
773  else
774  global_purity = static_cast<double>(global_nshared) / static_cast<double>(n_selected_hits);
775  } else
776  global_purity = 0;
777 
778  bool trackerOk = false;
779  if (n_tracker_selected_hits != 0) {
780  if (tracker_quality > tracker_quality_cut)
781  trackerOk = true;
782 
783  tracker_purity = static_cast<double>(tracker_nshared) / static_cast<double>(n_tracker_selected_hits);
785  tracker_purity = static_cast<double>(tracker_nshared);
786 
787  if ((!AbsoluteNumberOfHits_track) && tracker_purity <= PurityCut_track)
788  trackerOk = false;
789 
790  // if a track has just 3 hits in the Tracker we require that all 3
791  // hits are shared
792  if (ThreeHitTracksAreSpecial && n_tracker_selected_hits == 3 && tracker_nshared < 3)
793  trackerOk = false;
794  }
795 
796  bool muonOk = false;
797  if (n_muon_selected_hits != 0) {
798  if (muon_quality > muon_quality_cut)
799  muonOk = true;
800 
801  muon_purity = static_cast<double>(muon_nshared) / static_cast<double>(n_muon_selected_hits);
803  muon_purity = static_cast<double>(muon_nshared);
804 
805  if ((!AbsoluteNumberOfHits_muon) && muon_purity <= PurityCut_muon)
806  muonOk = false;
807  }
808 
809  // (matchOk) has to account for different track types (tracker-only,
810  // standalone muons, global muons)
811  bool matchOk = trackerOk || muonOk;
812 
813  // only for global muons: match both tracker and muon stub unless
814  // (acceptOneStubMatchings==true)
815  if (!acceptOneStubMatchings && n_tracker_selected_hits != 0 && n_muon_selected_hits != 0)
816  matchOk = trackerOk && muonOk;
817 
818  if (matchOk) {
819  outputCollection[tpindex].push_back(IndexMatch(tindex, global_quality));
820  any_trackingParticle_matched = true;
821 
822  edm::LogVerbatim("MuonAssociatorByHitsHelper")
823  << "*************************************************************"
824  "***********************************************************"
825  << "\n"
826  << "TrackingParticle " << tpindex << ", q = " << (*trpart).charge() << ", p = " << (*trpart).p()
827  << ", pT = " << (*trpart).pt() << ", eta = " << (*trpart).eta() << ", phi = " << (*trpart).phi() << "\n"
828  << " pdg code = " << (*trpart).pdgId() << ", made of " << (*trpart).numberOfHits()
829  << " PSimHits, recounted " << n_global_simhits << " PSimHits"
830  << " (tracker:" << n_tracker_recounted_simhits << "/muons:" << n_muon_simhits << ")"
831  << ", from " << (*trpart).g4Tracks().size() << " SimTrack:";
832  for (TrackingParticle::g4t_iterator g4T = (*trpart).g4Track_begin(); g4T != (*trpart).g4Track_end(); ++g4T) {
833  edm::LogVerbatim("MuonAssociatorByHitsHelper")
834  << " Id:" << (*g4T).trackId() << "/Evt:(" << (*g4T).eventId().event() << ","
835  << (*g4T).eventId().bunchCrossing() << ")";
836  }
837  edm::LogVerbatim("MuonAssociatorByHitsHelper")
838  << "\t selected " << n_global_selected_simhits << " PSimHits"
839  << " (tracker:" << n_tracker_selected_simhits << "/muons:" << n_muon_selected_simhits << ")"
840  << "\n\t **MATCHED** with quality = " << global_quality << " (tracker: " << tracker_quality
841  << " / muon: " << muon_quality << ")"
842  << "\n\t and purity = " << global_purity << " (tracker: " << tracker_purity
843  << " / muon: " << muon_purity << ")"
844  << "\n\t N shared hits = " << global_nshared << " (tracker: " << tracker_nshared
845  << " / muon: " << muon_nshared << ")"
846  << "\n"
847  << " to: reco::Track " << tindex << ZeroHitMuon << "\n\t"
848  << " made of " << n_selected_hits << " RecHits (tracker:" << n_tracker_valid
849  << "/muons:" << n_muon_selected_hits << ")";
850  } else {
851  // print something only if this TrackingParticle shares some hits with
852  // the current reco::Track
853  if (global_nshared != 0) {
854  if (printRtS)
855  edm::LogVerbatim("MuonAssociatorByHitsHelper")
856  << "*********************************************************"
857  "*********************************************************"
858  "******"
859  << "\n"
860  << "TrackingParticle " << tpindex << ", q = " << (*trpart).charge() << ", p = " << (*trpart).p()
861  << ", pT = " << (*trpart).pt() << ", eta = " << (*trpart).eta() << ", phi = " << (*trpart).phi()
862  << "\n"
863  << " pdg code = " << (*trpart).pdgId() << ", made of " << (*trpart).numberOfHits()
864  << " PSimHits, recounted " << n_global_simhits << " PSimHits"
865  << " (tracker:" << n_tracker_recounted_simhits << "/muons:" << n_muon_simhits << ")"
866  << ", from " << (*trpart).g4Tracks().size() << " SimTrack:";
867  for (TrackingParticle::g4t_iterator g4T = (*trpart).g4Track_begin(); g4T != (*trpart).g4Track_end();
868  ++g4T) {
869  if (printRtS)
870  edm::LogVerbatim("MuonAssociatorByHitsHelper")
871  << " Id:" << (*g4T).trackId() << "/Evt:(" << (*g4T).eventId().event() << ","
872  << (*g4T).eventId().bunchCrossing() << ")";
873  }
874  if (printRtS)
875  edm::LogVerbatim("MuonAssociatorByHitsHelper")
876  << "\t selected " << n_global_selected_simhits << " PSimHits"
877  << " (tracker:" << n_tracker_selected_simhits << "/muons:" << n_muon_selected_simhits << ")"
878  << "\n\t NOT matched to reco::Track " << tindex << ZeroHitMuon
879  << " with quality = " << global_quality << " (tracker: " << tracker_quality
880  << " / muon: " << muon_quality << ")"
881  << "\n\t and purity = " << global_purity << " (tracker: " << tracker_purity
882  << " / muon: " << muon_purity << ")"
883  << "\n\t N shared hits = " << global_nshared << " (tracker: " << tracker_nshared
884  << " / muon: " << muon_nshared << ")";
885  }
886  }
887  } // loop over TrackingParticle's
888  } // if(n_matching_simhits != 0)
889  } // loop over reco Tracks
890 
891  if (!any_trackingParticle_matched) {
892  edm::LogVerbatim("MuonAssociatorByHitsHelper")
893  << "\n"
894  << "*******************************************************************"
895  "*****************************************************"
896  << "\n NO TrackingParticle associated to ANY input reco::Track ! \n"
897  << "*******************************************************************"
898  "*****************************************************"
899  << "\n";
900  } else {
901  edm::LogVerbatim("MuonAssociatorByHitsHelper")
902  << "*******************************************************************"
903  "*****************************************************"
904  << "\n";
905  }
906 
907  for (IndexAssociation::iterator it = outputCollection.begin(), ed = outputCollection.end(); it != ed; ++it) {
908  std::sort(it->second.begin(), it->second.end());
909  }
910  return outputCollection;
911 }
912 
914  MapOfMatchedIds &muon_matchedIds_valid,
915  MapOfMatchedIds &tracker_matchedIds_INVALID,
916  MapOfMatchedIds &muon_matchedIds_INVALID,
917  int &n_tracker_valid,
918  int &n_dt_valid,
919  int &n_csc_valid,
920  int &n_rpc_valid,
921  int &n_gem_valid,
922  int &n_tracker_matched_valid,
923  int &n_dt_matched_valid,
924  int &n_csc_matched_valid,
925  int &n_rpc_matched_valid,
926  int &n_gem_matched_valid,
927  int &n_tracker_INVALID,
928  int &n_dt_INVALID,
929  int &n_csc_INVALID,
930  int &n_rpc_INVALID,
931  int &n_gem_INVALID,
932  int &n_tracker_matched_INVALID,
933  int &n_dt_matched_INVALID,
934  int &n_csc_matched_INVALID,
935  int &n_rpc_matched_INVALID,
936  int &n_gem_matched_INVALID,
939  const TrackerHitAssociator *trackertruth,
940  const DTHitAssociator &dttruth,
941  const CSCHitAssociator &csctruth,
942  const RPCHitAssociator &rpctruth,
943  const GEMHitAssociator &gemtruth,
944  bool printRtS,
945  const TrackerTopology *tTopo) const
946 
947 {
948  tracker_matchedIds_valid.clear();
949  muon_matchedIds_valid.clear();
950 
951  tracker_matchedIds_INVALID.clear();
952  muon_matchedIds_INVALID.clear();
953 
954  n_tracker_valid = 0;
955  n_dt_valid = 0;
956  n_csc_valid = 0;
957  n_rpc_valid = 0;
958  n_gem_valid = 0;
959 
960  n_tracker_matched_valid = 0;
961  n_dt_matched_valid = 0;
962  n_csc_matched_valid = 0;
963  n_rpc_matched_valid = 0;
964  n_gem_matched_valid = 0;
965 
966  n_tracker_INVALID = 0;
967  n_dt_INVALID = 0;
968  n_csc_INVALID = 0;
969  n_rpc_INVALID = 0;
970  n_gem_INVALID = 0;
971 
972  n_tracker_matched_INVALID = 0;
973  n_dt_matched_INVALID = 0;
974  n_csc_matched_INVALID = 0;
975  n_rpc_matched_INVALID = 0;
976  n_gem_matched_INVALID = 0;
977 
978  std::vector<SimHitIdpr> SimTrackIds;
979 
980  // main loop on TrackingRecHits
981  int iloop = 0;
982  int iH = -1;
983  for (trackingRecHit_iterator it = begin; it != end; it++, iloop++) {
984  stringstream hit_index;
985  hit_index << iloop;
986 
987  const TrackingRecHit *hitp = getHitPtr(it);
988  DetId geoid = hitp->geographicalId();
989 
990  unsigned int detid = geoid.rawId();
991  stringstream detector_id;
992  detector_id << detid;
993 
994  string hitlog = "TrackingRecHit " + hit_index.str();
995  string wireidlog;
996  std::vector<string> DTSimHits;
997 
998  DetId::Detector det = geoid.det();
999  int subdet = geoid.subdetId();
1000 
1001  bool valid_Hit = hitp->isValid();
1002 
1003  // Si-Tracker Hits
1004  if (det == DetId::Tracker && UseTracker) {
1005  stringstream detector_id;
1006  detector_id << tTopo->print(detid);
1007 
1008  if (valid_Hit)
1009  hitlog = hitlog + " -Tracker - detID = " + detector_id.str();
1010  else
1011  hitlog = hitlog + " *** INVALID ***" + " -Tracker - detID = " + detector_id.str();
1012 
1013  iH++;
1014  SimTrackIds = trackertruth->associateHitId(*hitp);
1015 
1016  if (valid_Hit) {
1017  n_tracker_valid++;
1018 
1019  if (!SimTrackIds.empty()) {
1020  n_tracker_matched_valid++;
1021  // tracker_matchedIds_valid[iH] = SimTrackIds;
1022  tracker_matchedIds_valid.push_back(new uint_SimHitIdpr_pair(iH, SimTrackIds));
1023  }
1024  } else {
1025  n_tracker_INVALID++;
1026 
1027  if (!SimTrackIds.empty()) {
1028  n_tracker_matched_INVALID++;
1029  // tracker_matchedIds_INVALID[iH] = SimTrackIds;
1030  tracker_matchedIds_INVALID.push_back(new uint_SimHitIdpr_pair(iH, SimTrackIds));
1031  }
1032  }
1033  }
1034  // Muon detector Hits
1035  else if (det == DetId::Muon && UseMuon) {
1036  // DT Hits
1037  if (subdet == MuonSubdetId::DT) {
1038  DTWireId dtdetid = DTWireId(detid);
1039  stringstream dt_detector_id;
1040  dt_detector_id << dtdetid;
1041  if (valid_Hit)
1042  hitlog = hitlog + " -Muon DT - detID = " + dt_detector_id.str();
1043  else
1044  hitlog = hitlog + " *** INVALID ***" + " -Muon DT - detID = " + dt_detector_id.str();
1045 
1046  const DTRecHit1D *dtrechit = dynamic_cast<const DTRecHit1D *>(hitp);
1047 
1048  // single DT hits
1049  if (dtrechit) {
1050  iH++;
1051  SimTrackIds = dttruth.associateDTHitId(dtrechit);
1052 
1053  if (valid_Hit) {
1054  n_dt_valid++;
1055 
1056  if (!SimTrackIds.empty()) {
1057  n_dt_matched_valid++;
1058  // muon_matchedIds_valid[iH] = SimTrackIds;
1059  muon_matchedIds_valid.push_back(new uint_SimHitIdpr_pair(iH, SimTrackIds));
1060  }
1061  } else {
1062  n_dt_INVALID++;
1063 
1064  if (!SimTrackIds.empty()) {
1065  n_dt_matched_INVALID++;
1066  // muon_matchedIds_INVALID[iH] = SimTrackIds;
1067  muon_matchedIds_INVALID.push_back(new uint_SimHitIdpr_pair(iH, SimTrackIds));
1068  }
1069  }
1070 
1071  if (dumpDT) {
1072  DTWireId wireid = dtrechit->wireId();
1073  stringstream wid;
1074  wid << wireid;
1075  std::vector<PSimHit> dtSimHits = dttruth.associateHit(*hitp);
1076 
1077  stringstream ndthits;
1078  ndthits << dtSimHits.size();
1079  wireidlog = "\t DTWireId :" + wid.str() + ", " + ndthits.str() + " associated PSimHit :";
1080 
1081  for (unsigned int j = 0; j < dtSimHits.size(); j++) {
1082  stringstream index;
1083  index << j;
1084  stringstream simhit;
1085  simhit << dtSimHits[j];
1086  string simhitlog = "\t\t PSimHit " + index.str() + ": " + simhit.str();
1087  DTSimHits.push_back(simhitlog);
1088  }
1089  } // if (dumpDT)
1090  }
1091 
1092  // DT segments
1093  else {
1094  const DTRecSegment4D *dtsegment = dynamic_cast<const DTRecSegment4D *>(hitp);
1095 
1096  if (dtsegment) {
1097  std::vector<const TrackingRecHit *> componentHits, phiHits, zHits;
1098  if (dtsegment->hasPhi()) {
1099  phiHits = dtsegment->phiSegment()->recHits();
1100  componentHits.insert(componentHits.end(), phiHits.begin(), phiHits.end());
1101  }
1102  if (dtsegment->hasZed()) {
1103  zHits = dtsegment->zSegment()->recHits();
1104  componentHits.insert(componentHits.end(), zHits.begin(), zHits.end());
1105  }
1106  if (printRtS)
1107  edm::LogVerbatim("MuonAssociatorByHitsHelper")
1108  << "\n\t this TrackingRecHit is a DTRecSegment4D with " << componentHits.size()
1109  << " hits (phi:" << phiHits.size() << ", z:" << zHits.size() << ")";
1110 
1111  SimTrackIds.clear();
1112  std::vector<SimHitIdpr> i_SimTrackIds;
1113  int i_compHit = 0;
1114  for (std::vector<const TrackingRecHit *>::const_iterator ithit = componentHits.begin();
1115  ithit != componentHits.end();
1116  ++ithit) {
1117  i_compHit++;
1118 
1119  const DTRecHit1D *dtrechit1D = dynamic_cast<const DTRecHit1D *>(*ithit);
1120 
1121  i_SimTrackIds.clear();
1122  if (dtrechit1D) {
1123  iH++;
1124  i_SimTrackIds = dttruth.associateDTHitId(dtrechit1D);
1125 
1126  if (valid_Hit) {
1127  // validity check is on the segment, but hits are counted
1128  // one-by-one
1129  n_dt_valid++;
1130 
1131  if (!i_SimTrackIds.empty()) {
1132  n_dt_matched_valid++;
1133  // muon_matchedIds_valid[iH] = i_SimTrackIds;
1134  muon_matchedIds_valid.push_back(new uint_SimHitIdpr_pair(iH, i_SimTrackIds));
1135  }
1136  } else {
1137  n_dt_INVALID++;
1138 
1139  if (!i_SimTrackIds.empty()) {
1140  n_dt_matched_INVALID++;
1141  // muon_matchedIds_INVALID[iH] = i_SimTrackIds;
1142  muon_matchedIds_INVALID.push_back(new uint_SimHitIdpr_pair(iH, i_SimTrackIds));
1143  }
1144  }
1145  } else if (printRtS)
1146  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "*** WARNING in "
1147  "MuonAssociatorByHitsHelper::getMatchedIds, null "
1148  "dynamic_cast of a DT TrackingRecHit !";
1149 
1150  unsigned int i_detid = (*ithit)->geographicalId().rawId();
1151  DTWireId i_dtdetid = DTWireId(i_detid);
1152 
1153  stringstream i_dt_detector_id;
1154  i_dt_detector_id << i_dtdetid;
1155 
1156  stringstream i_ss;
1157  i_ss << "\t\t hit " << i_compHit << " -Muon DT - detID = " << i_dt_detector_id.str();
1158 
1159  string i_hitlog = i_ss.str();
1160  i_hitlog = i_hitlog + write_matched_simtracks(i_SimTrackIds);
1161  if (printRtS)
1162  edm::LogVerbatim("MuonAssociatorByHitsHelper") << i_hitlog;
1163 
1164  SimTrackIds.insert(SimTrackIds.end(), i_SimTrackIds.begin(), i_SimTrackIds.end());
1165  }
1166  } // if (dtsegment)
1167 
1168  else if (printRtS)
1169  edm::LogVerbatim("MuonAssociatorByHitsHelper")
1170  << "*** WARNING in MuonAssociatorByHitsHelper::getMatchedIds, "
1171  "DT TrackingRecHit is neither DTRecHit1D nor DTRecSegment4D "
1172  "! ";
1173  }
1174  }
1175 
1176  // CSC Hits
1177  else if (subdet == MuonSubdetId::CSC) {
1178  CSCDetId cscdetid = CSCDetId(detid);
1179  stringstream csc_detector_id;
1180  csc_detector_id << cscdetid;
1181  if (valid_Hit)
1182  hitlog = hitlog + " -Muon CSC- detID = " + csc_detector_id.str();
1183  else
1184  hitlog = hitlog + " *** INVALID ***" + " -Muon CSC- detID = " + csc_detector_id.str();
1185 
1186  const CSCRecHit2D *cscrechit = dynamic_cast<const CSCRecHit2D *>(hitp);
1187 
1188  // single CSC hits
1189  if (cscrechit) {
1190  iH++;
1191  SimTrackIds = csctruth.associateCSCHitId(cscrechit);
1192 
1193  if (valid_Hit) {
1194  n_csc_valid++;
1195 
1196  if (!SimTrackIds.empty()) {
1197  n_csc_matched_valid++;
1198  // muon_matchedIds_valid[iH] = SimTrackIds;
1199  muon_matchedIds_valid.push_back(new uint_SimHitIdpr_pair(iH, SimTrackIds));
1200  }
1201  } else {
1202  n_csc_INVALID++;
1203 
1204  if (!SimTrackIds.empty()) {
1205  n_csc_matched_INVALID++;
1206  // muon_matchedIds_INVALID[iH] = SimTrackIds;
1207  muon_matchedIds_INVALID.push_back(new uint_SimHitIdpr_pair(iH, SimTrackIds));
1208  }
1209  }
1210  }
1211 
1212  // CSC segments
1213  else {
1214  const CSCSegment *cscsegment = dynamic_cast<const CSCSegment *>(hitp);
1215 
1216  if (cscsegment) {
1217  std::vector<const TrackingRecHit *> componentHits = cscsegment->recHits();
1218  if (printRtS)
1219  edm::LogVerbatim("MuonAssociatorByHitsHelper")
1220  << "\n\t this TrackingRecHit is a CSCSegment with " << componentHits.size() << " hits";
1221 
1222  SimTrackIds.clear();
1223  std::vector<SimHitIdpr> i_SimTrackIds;
1224  int i_compHit = 0;
1225  for (std::vector<const TrackingRecHit *>::const_iterator ithit = componentHits.begin();
1226  ithit != componentHits.end();
1227  ++ithit) {
1228  i_compHit++;
1229 
1230  const CSCRecHit2D *cscrechit2D = dynamic_cast<const CSCRecHit2D *>(*ithit);
1231 
1232  i_SimTrackIds.clear();
1233  if (cscrechit2D) {
1234  iH++;
1235  i_SimTrackIds = csctruth.associateCSCHitId(cscrechit2D);
1236 
1237  if (valid_Hit) {
1238  // validity check is on the segment, but hits are counted
1239  // one-by-one
1240  n_csc_valid++;
1241 
1242  if (!i_SimTrackIds.empty()) {
1243  n_csc_matched_valid++;
1244  // muon_matchedIds_valid[iH] = i_SimTrackIds;
1245  muon_matchedIds_valid.push_back(new uint_SimHitIdpr_pair(iH, i_SimTrackIds));
1246  }
1247  } else {
1248  n_csc_INVALID++;
1249 
1250  if (!i_SimTrackIds.empty()) {
1251  n_csc_matched_INVALID++;
1252  // muon_matchedIds_INVALID[iH] = i_SimTrackIds;
1253  muon_matchedIds_INVALID.push_back(new uint_SimHitIdpr_pair(iH, i_SimTrackIds));
1254  }
1255  }
1256  } else if (printRtS)
1257  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "*** WARNING in "
1258  "MuonAssociatorByHitsHelper::getMatchedIds, null "
1259  "dynamic_cast of a CSC TrackingRecHit !";
1260 
1261  unsigned int i_detid = (*ithit)->geographicalId().rawId();
1262  CSCDetId i_cscdetid = CSCDetId(i_detid);
1263 
1264  stringstream i_csc_detector_id;
1265  i_csc_detector_id << i_cscdetid;
1266 
1267  stringstream i_ss;
1268  i_ss << "\t\t hit " << i_compHit << " -Muon CSC- detID = " << i_csc_detector_id.str();
1269 
1270  string i_hitlog = i_ss.str();
1271  i_hitlog = i_hitlog + write_matched_simtracks(i_SimTrackIds);
1272  if (printRtS)
1273  edm::LogVerbatim("MuonAssociatorByHitsHelper") << i_hitlog;
1274 
1275  SimTrackIds.insert(SimTrackIds.end(), i_SimTrackIds.begin(), i_SimTrackIds.end());
1276  }
1277  } // if (cscsegment)
1278 
1279  else if (printRtS)
1280  edm::LogVerbatim("MuonAssociatorByHitsHelper")
1281  << "*** WARNING in MuonAssociatorByHitsHelper::getMatchedIds, "
1282  "CSC TrackingRecHit is neither CSCRecHit2D nor CSCSegment "
1283  "! ";
1284  }
1285  }
1286 
1287  // RPC Hits
1288  else if (subdet == MuonSubdetId::RPC) {
1289  RPCDetId rpcdetid = RPCDetId(detid);
1290  stringstream rpc_detector_id;
1291  rpc_detector_id << rpcdetid;
1292  if (valid_Hit)
1293  hitlog = hitlog + " -Muon RPC- detID = " + rpc_detector_id.str();
1294  else
1295  hitlog = hitlog + " *** INVALID ***" + " -Muon RPC- detID = " + rpc_detector_id.str();
1296 
1297  iH++;
1298  SimTrackIds = rpctruth.associateRecHit(*hitp);
1299 
1300  if (valid_Hit) {
1301  n_rpc_valid++;
1302 
1303  if (!SimTrackIds.empty()) {
1304  n_rpc_matched_valid++;
1305  // muon_matchedIds_valid[iH] = SimTrackIds;
1306  muon_matchedIds_valid.push_back(new uint_SimHitIdpr_pair(iH, SimTrackIds));
1307  }
1308  } else {
1309  n_rpc_INVALID++;
1310 
1311  if (!SimTrackIds.empty()) {
1312  n_rpc_matched_INVALID++;
1313  // muon_matchedIds_INVALID[iH] = SimTrackIds;
1314  muon_matchedIds_INVALID.push_back(new uint_SimHitIdpr_pair(iH, SimTrackIds));
1315  }
1316  }
1317  }
1318 
1319  // GEM Hits
1320  else if (subdet == MuonSubdetId::GEM) {
1321  GEMDetId gemdetid = GEMDetId(detid);
1322  stringstream gem_detector_id;
1323  gem_detector_id << gemdetid;
1324  if (valid_Hit)
1325  hitlog = hitlog + " -Muon GEM- detID = " + gem_detector_id.str();
1326  else
1327  hitlog = hitlog + " *** INVALID ***" + " -Muon GEM- detID = " + gem_detector_id.str();
1328 
1329  const GEMRecHit *gemrechit = dynamic_cast<const GEMRecHit *>(hitp);
1330  if (gemrechit) {
1331  iH++;
1332  SimTrackIds = gemtruth.associateRecHit(gemrechit);
1333 
1334  if (valid_Hit) {
1335  n_gem_valid++;
1336 
1337  if (!SimTrackIds.empty()) {
1338  n_gem_matched_valid++;
1339  // muon_matchedIds_valid[iH] = SimTrackIds;
1340  muon_matchedIds_valid.push_back(new uint_SimHitIdpr_pair(iH, SimTrackIds));
1341  }
1342  } else {
1343  n_gem_INVALID++;
1344 
1345  if (!SimTrackIds.empty()) {
1346  n_gem_matched_INVALID++;
1347  // muon_matchedIds_INVALID[iH] = SimTrackIds;
1348  muon_matchedIds_INVALID.push_back(new uint_SimHitIdpr_pair(iH, SimTrackIds));
1349  }
1350  }
1351  } else {
1352  const GEMSegment *gemsegment = dynamic_cast<const GEMSegment *>(hitp);
1353  if (gemsegment) {
1354  std::vector<const TrackingRecHit *> componentHits = gemsegment->recHits();
1355  if (printRtS)
1356  edm::LogVerbatim("MuonAssociatorByHitsHelper")
1357  << "\n\t this TrackingRecHit is a GEMSegment with " << componentHits.size() << " hits";
1358 
1359  SimTrackIds.clear();
1360  std::vector<SimHitIdpr> i_SimTrackIds;
1361  int i_compHit = 0;
1362 
1363  for (auto const &ithit : componentHits) {
1364  i_compHit++;
1365 
1366  const GEMRecHit *gemrechitseg = dynamic_cast<const GEMRecHit *>(ithit);
1367 
1368  i_SimTrackIds.clear();
1369  if (gemrechitseg) {
1370  iH++;
1371  i_SimTrackIds = gemtruth.associateRecHit(gemrechitseg);
1372 
1373  if (valid_Hit) {
1374  // validity check is on the segment, but hits are counted
1375  // one-by-one
1376  n_gem_valid++;
1377 
1378  if (!i_SimTrackIds.empty()) {
1379  n_gem_matched_valid++;
1380  // muon_matchedIds_valid[iH] = i_SimTrackIds;
1381  muon_matchedIds_valid.push_back(new uint_SimHitIdpr_pair(iH, i_SimTrackIds));
1382  }
1383  } else {
1384  n_gem_INVALID++;
1385 
1386  if (!i_SimTrackIds.empty()) {
1387  n_gem_matched_INVALID++;
1388  // muon_matchedIds_INVALID[iH] = i_SimTrackIds;
1389  muon_matchedIds_INVALID.push_back(new uint_SimHitIdpr_pair(iH, i_SimTrackIds));
1390  }
1391  }
1392  } else if (printRtS)
1393  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "*** WARNING in "
1394  "MuonAssociatorByHitsHelper::getMatchedIds, null "
1395  "dynamic_cast of a GEM TrackingRecHit !";
1396 
1397  if (printRtS) {
1398  unsigned int i_detid = ithit->geographicalId().rawId();
1399  GEMDetId i_gemdetid = GEMDetId(i_detid);
1400 
1401  string i_hitlog = std::to_string(i_gemdetid);
1402  i_hitlog = i_hitlog + write_matched_simtracks(i_SimTrackIds);
1403  edm::LogVerbatim("MuonAssociatorByHitsHelper") << i_hitlog;
1404  }
1405 
1406  SimTrackIds.insert(SimTrackIds.end(), i_SimTrackIds.begin(), i_SimTrackIds.end());
1407  }
1408  } // if (gemsegment)
1409 
1410  } // if (gemrechit
1411  } else if (printRtS)
1412  edm::LogVerbatim("MuonAssociatorByHitsHelper")
1413  << "TrackingRecHit " << iloop << " *** WARNING *** Unexpected Hit from Detector = " << det;
1414  } // end if (det == DetId::Muon && UseMuon)
1415  else
1416  continue;
1417 
1418  hitlog = hitlog + write_matched_simtracks(SimTrackIds);
1419 
1420  if (printRtS)
1421  edm::LogVerbatim("MuonAssociatorByHitsHelper") << hitlog;
1422  if (printRtS && dumpDT && det == DetId::Muon && subdet == MuonSubdetId::DT) {
1423  edm::LogVerbatim("MuonAssociatorByHitsHelper") << wireidlog;
1424  for (unsigned int j = 0; j < DTSimHits.size(); j++) {
1425  edm::LogVerbatim("MuonAssociatorByHitsHelper") << DTSimHits[j];
1426  }
1427  }
1428 
1429  } // trackingRecHit loop
1430 }
1431 
1433  TrackingParticleCollection::const_iterator trpart) const {
1434  int nshared = 0;
1435  const std::vector<SimTrack> &g4Tracks = trpart->g4Tracks();
1436 
1437  // map is indexed over the rechits of the reco::Track (no double-countings
1438  // allowed)
1439  for (MapOfMatchedIds::const_iterator iRecH = matchedIds.begin(); iRecH != matchedIds.end(); ++iRecH) {
1440  // vector of associated simhits associated to the current rechit
1441  std::vector<SimHitIdpr> const &SimTrackIds = (*iRecH).second;
1442 
1443  bool found = false;
1444 
1445  for (const auto &iSimH : SimTrackIds) {
1446  uint32_t simtrackId = iSimH.first;
1447  EncodedEventId evtId = iSimH.second;
1448 
1449  // look for shared hits with the given TrackingParticle (looping over
1450  // component SimTracks)
1451  for (const auto &simtrack : g4Tracks) {
1452  if (simtrack.trackId() == simtrackId && simtrack.eventId() == evtId) {
1453  found = true;
1454  break;
1455  }
1456  }
1457 
1458  if (found) {
1459  nshared++;
1460  break;
1461  }
1462  }
1463  }
1464 
1465  return nshared;
1466 }
1467 
1468 std::string MuonAssociatorByHitsHelper::write_matched_simtracks(const std::vector<SimHitIdpr> &SimTrackIds) const {
1469  if (SimTrackIds.empty())
1470  return " *** UNMATCHED ***";
1471 
1472  string hitlog(" matched to SimTrack");
1473 
1474  for (size_t j = 0; j < SimTrackIds.size(); j++) {
1475  char buf[64];
1476  snprintf(buf,
1477  64,
1478  " Id:%i/Evt:(%i,%i) ",
1479  SimTrackIds[j].first,
1480  SimTrackIds[j].second.event(),
1481  SimTrackIds[j].second.bunchCrossing());
1482  hitlog += buf;
1483  }
1484  return hitlog;
1485 }
static constexpr int GEM
Definition: MuonSubdetId.h:15
std::vector< TrackingParticle > TrackingParticleCollection
std::string dump(unsigned int indent=0) const
std::function< void(const TrackHitsCollection &, const TrackingParticleCollection &)> diagnostics_
std::string write_matched_simtracks(const std::vector< SimHitIdpr > &) const
const DTChamberRecSegment2D * phiSegment() const
The superPhi segment: 0 if no phi projection available.
std::map< size_t, std::vector< IndexMatch > > IndexAssociation
int getShared(MapOfMatchedIds &matchedIds, TrackingParticleCollection::const_iterator trpart) const
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:50
std::vector< SimHitIdpr > associateCSCHitId(const CSCRecHit2D *) const
std::pair< unsigned int, std::vector< SimHitIdpr > > uint_SimHitIdpr_pair
std::string print(DetId detid) const
std::vector< PSimHit > associateHit(const TrackingRecHit &hit) const
IndexAssociation associateSimToRecoIndices(const TrackHitsCollection &, const edm::RefVector< TrackingParticleCollection > &, Resources const &) const
U second(std::pair< T, U > const &p)
std::vector< SimHitIdpr > associateRecHit(const TrackingRecHit &hit) const
std::vector< std::pair< trackingRecHit_iterator, trackingRecHit_iterator > > TrackHitsCollection
MuonAssociatorByHitsHelper(const edm::ParameterSet &conf)
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:41
#define end
Definition: vmac.h:39
bool hasPhi() const
Does it have the Phi projection?
std::vector< SimTrack >::const_iterator g4t_iterator
bool hasZed() const
Does it have the Z projection?
Definition: DetId.h:18
const DTSLRecSegment2D * zSegment() const
The Z segment: 0 if not zed projection available.
std::vector< const TrackingRecHit * > recHits() const override
Access to component RecHits (if any)
Definition: CSCSegment.cc:30
std::vector< SimHitIdpr > associateRecHit(const GEMRecHit *gemrechit) const
Detector
Definition: DetId.h:26
static constexpr int RPC
Definition: MuonSubdetId.h:14
bool isValid() const
IndexAssociation associateRecoToSimIndices(const TrackHitsCollection &, const edm::RefVector< TrackingParticleCollection > &, Resources const &) const
std::vector< SimHitIdpr > associateHitId(const TrackingRecHit &thit) const
fixed size matrix
#define begin
Definition: vmac.h:32
boost::ptr_vector< uint_SimHitIdpr_pair > MapOfMatchedIds
const TrackingRecHit * getHitPtr(edm::OwnVector< TrackingRecHit >::const_iterator iter) const
size_type size() const
Size of the RefVector.
Definition: RefVector.h:107
std::vector< const TrackingRecHit * > recHits() const override
Access to component RecHits (if any)
DetId geographicalId() const
static constexpr int DT
Definition: MuonSubdetId.h:12
std::vector< const TrackingRecHit * > recHits() const override
Access to component RecHits (if any)
Definition: GEMSegment.cc:72
static constexpr int CSC
Definition: MuonSubdetId.h:13
std::vector< SimHitIdpr > associateDTHitId(const DTRecHit1D *dtrechit) const
void getMatchedIds(MapOfMatchedIds &tracker_matchedIds_valid, MapOfMatchedIds &muon_matchedIds_valid, MapOfMatchedIds &tracker_matchedIds_INVALID, MapOfMatchedIds &muon_matchedIds_INVALID, int &n_tracker_valid, int &n_dt_valid, int &n_csc_valid, int &n_rpc_valid, int &n_gem_valid, int &n_tracker_matched_valid, int &n_dt_matched_valid, int &n_csc_matched_valid, int &n_rpc_matched_valid, int &n_gem_matched_valid, int &n_tracker_INVALID, int &n_dt_INVALID, int &n_csc_INVALID, int &n_rpc_INVALID, int &n_gem_INVALID, int &n_tracker_matched_INVALID, int &n_dt_matched_INVALID, int &n_csc_matched_INVALID, int &n_rpc_matched_INVALID, int &n_gem_matched_INVALID, trackingRecHit_iterator begin, trackingRecHit_iterator end, const TrackerHitAssociator *trackertruth, const DTHitAssociator &dttruth, const CSCHitAssociator &csctruth, const RPCHitAssociator &rpctruth, const GEMHitAssociator &gemtruth, bool printRts, const TrackerTopology *) const
constexpr Detector det() const
get the detector field from this detid
Definition: DetId.h:39
DTWireId wireId() const
Return the wireId.
Definition: DTRecHit1D.h:107