CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
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  rejectBadGlobal(conf.getParameter<bool>("rejectBadGlobal")),
21  UseTracker(conf.getParameter<bool>("UseTracker")),
22  UseMuon(conf.getParameter<bool>("UseMuon")),
23  AbsoluteNumberOfHits_track(conf.getParameter<bool>("AbsoluteNumberOfHits_track")),
24  NHitCut_track(conf.getParameter<unsigned int>("NHitCut_track")),
25  EfficiencyCut_track(conf.getParameter<double>("EfficiencyCut_track")),
26  PurityCut_track(conf.getParameter<double>("PurityCut_track")),
27  AbsoluteNumberOfHits_muon(conf.getParameter<bool>("AbsoluteNumberOfHits_muon")),
28  NHitCut_muon(conf.getParameter<unsigned int>("NHitCut_muon")),
29  EfficiencyCut_muon(conf.getParameter<double>("EfficiencyCut_muon")),
30  PurityCut_muon(conf.getParameter<double>("PurityCut_muon")),
31  UsePixels(conf.getParameter<bool>("UsePixels")),
32  UseGrouped(conf.getParameter<bool>("UseGrouped")),
33  UseSplitting(conf.getParameter<bool>("UseSplitting")),
34  ThreeHitTracksAreSpecial(conf.getParameter<bool>("ThreeHitTracksAreSpecial")),
35  dumpDT(conf.getParameter<bool>("dumpDT")) {
36  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "\n constructing MuonAssociatorByHitsHelper" << conf.dump();
37 
38  // up to the user in the other cases - print a message
39  if (UseTracker)
40  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "\n UseTracker = TRUE : Tracker SimHits and RecHits WILL be "
41  "counted";
42  else
43  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "\n UseTracker = FALSE : Tracker SimHits and RecHits WILL NOT be "
44  "counted";
45 
46  // up to the user in the other cases - print a message
47  if (UseMuon)
48  edm::LogVerbatim("MuonAssociatorByHitsHelper") << " UseMuon = TRUE : Muon SimHits and RecHits WILL be counted";
49  else
50  edm::LogVerbatim("MuonAssociatorByHitsHelper")
51  << " UseMuon = FALSE : Muon SimHits and RecHits WILL NOT be counted" << endl;
52 
53  // check consistency of the configuration when allowing zero-hit muon matching
54  // (counting invalid hits)
55  if (includeZeroHitMuons) {
56  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "\n includeZeroHitMuons = TRUE"
57  << "\n ==> (re)set NHitCut_muon = 0, PurityCut_muon = 0, "
58  "EfficiencyCut_muon = 0"
59  << endl;
60  NHitCut_muon = 0;
61  PurityCut_muon = 0.;
62  EfficiencyCut_muon = 0.;
63  }
64 }
65 
67  const TrackHitsCollection &tC,
68  const edm::RefVector<TrackingParticleCollection> &TPCollectionH,
69  const Resources &resources) const {
70  auto tTopo = resources.tTopo_;
71  auto trackertruth = resources.trackerHitAssoc_;
72  auto const &csctruth = *resources.cscHitAssoc_;
73  auto const &dttruth = *resources.dtHitAssoc_;
74  auto const &rpctruth = *resources.rpcHitAssoc_;
75  auto const &gemtruth = *resources.gemHitAssoc_;
76 
77  int tracker_nshared = 0;
78  int muon_nshared = 0;
79  int global_nshared = 0;
80 
81  double tracker_quality = 0;
82  double tracker_quality_cut;
84  tracker_quality_cut = static_cast<double>(NHitCut_track);
85  else
86  tracker_quality_cut = PurityCut_track;
87 
88  double muon_quality = 0;
89  double muon_quality_cut;
91  muon_quality_cut = static_cast<double>(NHitCut_muon);
92  else
93  muon_quality_cut = PurityCut_muon;
94 
95  double global_quality = 0;
96 
97  MapOfMatchedIds tracker_matchedIds_valid, muon_matchedIds_valid;
98  MapOfMatchedIds tracker_matchedIds_INVALID, muon_matchedIds_INVALID;
99 
100  IndexAssociation outputCollection;
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  bool 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 tracks: match both tracker and muon stub (if acceptOneStubMatchings==false)
342  // depending on the muon selection reject tracks with only one stub (if rejectBadGlobal==true)
343  //
345  ((n_tracker_selected_hits != 0 && n_muon_selected_hits != 0) || rejectBadGlobal))
346  matchOk = trackerOk && muonOk;
347 
348  if (matchOk) {
349  outputCollection[tindex].push_back(IndexMatch(tpindex, global_quality));
350  this_track_matched = true;
351 
352  edm::LogVerbatim("MuonAssociatorByHitsHelper")
353  << "\n\t"
354  << " **MATCHED** with quality = " << global_quality << " (tracker: " << tracker_quality
355  << " / muon: " << muon_quality << ")"
356  << "\n\t"
357  << " N shared hits = " << global_nshared << " (tracker: " << tracker_nshared
358  << " / muon: " << muon_nshared << ")"
359  << "\n"
360  << " to TrackingParticle: q = " << (*trpart).charge() << ", p = " << (*trpart).p()
361  << ", pT = " << (*trpart).pt() << ", eta = " << (*trpart).eta() << ", phi = " << (*trpart).phi() << "\n\t"
362  << " pdg code = " << (*trpart).pdgId() << ", made of " << (*trpart).numberOfHits() << " PSimHits"
363  << " from " << (*trpart).g4Tracks().size() << " SimTrack:";
364  for (TrackingParticle::g4t_iterator g4T = (*trpart).g4Track_begin(); g4T != (*trpart).g4Track_end(); ++g4T) {
365  edm::LogVerbatim("MuonAssociatorByHitsHelper")
366  << "\t"
367  << " Id:" << (*g4T).trackId() << "/Evt:(" << (*g4T).eventId().event() << ","
368  << (*g4T).eventId().bunchCrossing() << ")";
369  }
370  } else {
371  // print something only if this TrackingParticle shares some hits with
372  // the current reco::Track
373  if (global_nshared != 0)
374  LogTrace("MuonAssociatorByHitsHelper")
375  << "\n\t"
376  << " NOT matched to TrackingParticle, with quality = " << global_quality
377  << " (tracker: " << tracker_quality << " / muon: " << muon_quality << ")"
378  << "\n"
379  << " N shared hits = " << global_nshared << " (tracker: " << tracker_nshared
380  << " / muon: " << muon_nshared << ")";
381  }
382 
383  } // loop over TrackingParticle
384 
385  if (!this_track_matched) {
386  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "\n"
387  << " NOT matched to any TrackingParticle";
388  }
389 
390  edm::LogVerbatim("MuonAssociatorByHitsHelper")
391  << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
392  "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
393  << "\n";
394 
395  } // if(n_matching_simhits != 0)
396 
397  } // loop over reco::Track
398 
399  if (tC.empty())
400  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "0 reconstructed tracks (-->> 0 associated !)";
401 
402  for (IndexAssociation::iterator it = outputCollection.begin(), ed = outputCollection.end(); it != ed; ++it) {
403  std::sort(it->second.begin(), it->second.end());
404  }
405  return outputCollection;
406 }
407 
409  const TrackHitsCollection &tC,
410  const edm::RefVector<TrackingParticleCollection> &TPCollectionH,
411  const Resources &resources) const {
412  auto tTopo = resources.tTopo_;
413  auto trackertruth = resources.trackerHitAssoc_;
414  auto const &csctruth = *resources.cscHitAssoc_;
415  auto const &dttruth = *resources.dtHitAssoc_;
416  auto const &rpctruth = *resources.rpcHitAssoc_;
417  auto const &gemtruth = *resources.gemHitAssoc_;
418 
419  int tracker_nshared = 0;
420  int muon_nshared = 0;
421  int global_nshared = 0;
422 
423  double tracker_quality = 0;
424  double tracker_quality_cut;
426  tracker_quality_cut = static_cast<double>(NHitCut_track);
427  else
428  tracker_quality_cut = EfficiencyCut_track;
429 
430  double muon_quality = 0;
431  double muon_quality_cut;
433  muon_quality_cut = static_cast<double>(NHitCut_muon);
434  else
435  muon_quality_cut = EfficiencyCut_muon;
436 
437  double global_quality = 0;
438 
439  double tracker_purity = 0;
440  double muon_purity = 0;
441  double global_purity = 0;
442 
443  MapOfMatchedIds tracker_matchedIds_valid, muon_matchedIds_valid;
444  MapOfMatchedIds tracker_matchedIds_INVALID, muon_matchedIds_INVALID;
445 
446  IndexAssociation outputCollection;
447 
448  bool printRtS = false;
449 
451  tPC.reserve(TPCollectionH.size());
452  for (auto const &ref : TPCollectionH) {
453  tPC.push_back(*ref);
454  }
455 
456  bool any_trackingParticle_matched = false;
457 
458  int tindex = 0;
459  for (TrackHitsCollection::const_iterator track = tC.begin(); track != tC.end(); track++, tindex++) {
460  if (printRtS)
461  edm::LogVerbatim("MuonAssociatorByHitsHelper")
462  << "\n"
463  << "reco::Track " << tindex << ", number of RecHits = " << (track->second - track->first) << "\n";
464 
465  tracker_matchedIds_valid.clear();
466  muon_matchedIds_valid.clear();
467 
468  tracker_matchedIds_INVALID.clear();
469  muon_matchedIds_INVALID.clear();
470 
471  int n_matching_simhits = 0;
472 
473  // all hits = valid +INVALID
474  int n_all = 0;
475  int n_tracker_all = 0;
476  int n_dt_all = 0;
477  int n_csc_all = 0;
478  int n_rpc_all = 0;
479  int n_gem_all = 0;
480 
481  int n_valid = 0;
482  int n_tracker_valid = 0;
483  int n_muon_valid = 0;
484  int n_dt_valid = 0;
485  int n_csc_valid = 0;
486  int n_rpc_valid = 0;
487  int n_gem_valid = 0;
488 
489  int n_tracker_matched_valid = 0;
490  int n_muon_matched_valid = 0;
491  int n_dt_matched_valid = 0;
492  int n_csc_matched_valid = 0;
493  int n_rpc_matched_valid = 0;
494  int n_gem_matched_valid = 0;
495 
496  int n_INVALID = 0;
497  int n_tracker_INVALID = 0;
498  int n_muon_INVALID = 0;
499  int n_dt_INVALID = 0;
500  int n_csc_INVALID = 0;
501  int n_rpc_INVALID = 0;
502  int n_gem_INVALID = 0;
503 
504  int n_tracker_matched_INVALID = 0;
505  int n_muon_matched_INVALID = 0;
506  int n_dt_matched_INVALID = 0;
507  int n_csc_matched_INVALID = 0;
508  int n_rpc_matched_INVALID = 0;
509  int n_gem_matched_INVALID = 0;
510 
511  getMatchedIds(tracker_matchedIds_valid,
512  muon_matchedIds_valid,
513  tracker_matchedIds_INVALID,
514  muon_matchedIds_INVALID,
515  n_tracker_valid,
516  n_dt_valid,
517  n_csc_valid,
518  n_rpc_valid,
519  n_gem_valid,
520  n_tracker_matched_valid,
521  n_dt_matched_valid,
522  n_csc_matched_valid,
523  n_rpc_matched_valid,
524  n_gem_matched_valid,
525  n_tracker_INVALID,
526  n_dt_INVALID,
527  n_csc_INVALID,
528  n_rpc_INVALID,
529  n_gem_INVALID,
530  n_tracker_matched_INVALID,
531  n_dt_matched_INVALID,
532  n_csc_matched_INVALID,
533  n_rpc_matched_INVALID,
534  n_gem_matched_INVALID,
535  track->first,
536  track->second,
537  trackertruth,
538  dttruth,
539  csctruth,
540  rpctruth,
541  gemtruth,
542  printRtS,
543  tTopo);
544 
545  n_matching_simhits = tracker_matchedIds_valid.size() + muon_matchedIds_valid.size() +
546  tracker_matchedIds_INVALID.size() + muon_matchedIds_INVALID.size();
547 
548  n_muon_valid = n_dt_valid + n_csc_valid + n_rpc_valid + n_gem_valid;
549  n_valid = n_tracker_valid + n_muon_valid;
550  n_muon_INVALID = n_dt_INVALID + n_csc_INVALID + n_rpc_INVALID + n_gem_INVALID;
551  n_INVALID = n_tracker_INVALID + n_muon_INVALID;
552 
553  // all used hits (valid+INVALID), defined by UseTracker, UseMuon
554  n_tracker_all = n_tracker_valid + n_tracker_INVALID;
555  n_dt_all = n_dt_valid + n_dt_INVALID;
556  n_csc_all = n_csc_valid + n_csc_INVALID;
557  n_rpc_all = n_rpc_valid + n_rpc_INVALID;
558  n_gem_all = n_gem_valid + n_gem_INVALID;
559  n_all = n_valid + n_INVALID;
560 
561  n_muon_matched_valid = n_dt_matched_valid + n_csc_matched_valid + n_rpc_matched_valid + n_gem_matched_valid;
562  n_muon_matched_INVALID =
563  n_dt_matched_INVALID + n_csc_matched_INVALID + n_rpc_matched_INVALID + n_gem_matched_INVALID;
564 
565  // selected hits are set initially to valid hits
566  int n_tracker_selected_hits = n_tracker_valid;
567  int n_muon_selected_hits = n_muon_valid;
568  int n_dt_selected_hits = n_dt_valid;
569  int n_csc_selected_hits = n_csc_valid;
570  int n_rpc_selected_hits = n_rpc_valid;
571  int n_gem_selected_hits = n_gem_valid;
572 
573  // matched hits are a subsample of the selected hits
574  int n_tracker_matched = n_tracker_matched_valid;
575  int n_muon_matched = n_muon_matched_valid;
576  int n_dt_matched = n_dt_matched_valid;
577  int n_csc_matched = n_csc_matched_valid;
578  int n_rpc_matched = n_rpc_matched_valid;
579  int n_gem_matched = n_gem_matched_valid;
580 
581  std::string InvMuonHits, ZeroHitMuon;
582 
583  if (includeZeroHitMuons && n_muon_valid == 0 && n_muon_INVALID != 0) {
584  // selected muon hits = INVALID when (useZeroHitMuons == True) and track
585  // has no valid muon hits
586 
587  InvMuonHits = " ***INVALID MUON HITS***";
588  ZeroHitMuon = " ***ZERO-HIT MUON***";
589 
590  n_muon_selected_hits = n_muon_INVALID;
591  n_dt_selected_hits = n_dt_INVALID;
592  n_csc_selected_hits = n_csc_INVALID;
593  n_rpc_selected_hits = n_rpc_INVALID;
594  n_gem_selected_hits = n_gem_INVALID;
595 
596  n_muon_matched = n_muon_matched_INVALID;
597  n_dt_matched = n_dt_matched_INVALID;
598  n_csc_matched = n_csc_matched_INVALID;
599  n_rpc_matched = n_rpc_matched_INVALID;
600  n_gem_matched = n_gem_matched_INVALID;
601  }
602 
603  int n_selected_hits = n_tracker_selected_hits + n_muon_selected_hits;
604  int n_matched = n_tracker_matched + n_muon_matched;
605 
606  if (printRtS)
607  edm::LogVerbatim("MuonAssociatorByHitsHelper")
608  << "\n"
609  << "# TrackingRecHits: " << (track->second - track->first) << "\n"
610  << "# used RecHits = " << n_all << " (" << n_tracker_all << "/" << n_dt_all << "/" << n_csc_all << "/"
611  << n_rpc_all << "/" << n_gem_all << " in Tracker/DT/CSC/RPC/GEM)"
612  << ", obtained from " << n_matching_simhits << " SimHits"
613  << "\n"
614  << "# selected RecHits = " << n_selected_hits << " (" << n_tracker_selected_hits << "/" << n_dt_selected_hits
615  << "/" << n_csc_selected_hits << "/" << n_rpc_selected_hits << "/" << n_gem_selected_hits
616  << " in Tracker/DT/CSC/RPC/GEM)" << InvMuonHits << "\n"
617  << "# matched RecHits = " << n_matched << " (" << n_tracker_matched << "/" << n_dt_matched << "/"
618  << n_csc_matched << "/" << n_rpc_matched << "/" << n_gem_matched << " in Tracker/DT/CSC/RPC/GEM)";
619 
620  if (printRtS && n_all > 0 && n_matching_simhits == 0)
621  edm::LogVerbatim("MuonAssociatorByHitsHelper")
622  << "*** WARNING in MuonAssociatorByHitsHelper::associateSimToReco: "
623  "no matching PSimHit found for this reco::Track !";
624 
625  if (n_matching_simhits != 0) {
626  int tpindex = 0;
627  for (TrackingParticleCollection::const_iterator trpart = tPC.begin(); trpart != tPC.end(); ++trpart, ++tpindex) {
628  // int n_tracker_simhits = 0;
629  int n_tracker_recounted_simhits = 0;
630  int n_muon_simhits = 0;
631  int n_global_simhits = 0;
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  // adapt to new TP interface: this gives the total number of hits in
650  // tracker
651  // should reproduce the behaviour of UseGrouped=UseSplitting=.true.
652  n_tracker_recounted_simhits = trpart->numberOfTrackerHits();
653  // numberOfHits() gives the total number of hits (tracker + muons)
654  n_muon_simhits = trpart->numberOfHits() - trpart->numberOfTrackerHits();
655 
656  // Handle the case of TrackingParticles that don't have PSimHits inside,
657  // e.g. because they were made on RECOSIM only.
658  if (trpart->numberOfHits() == 0) {
659  // FIXME this can be made better, counting the digiSimLinks associated
660  // to this TP, but perhaps it's not worth it
661  // n_tracker_recounted_simhits = tracker_nshared;
662  // n_muon_simhits = muon_nshared;
663  // ---> on RECOSIM when the AbsoluteNumberOfHits_muon=True this always
664  // obtains quality=1,
665  // hence no sorting is possible in case of duplicate matchings
666  // ---> reset these variables to 1 so to keep the number of shared
667  // hits as ranking criterion
668  n_tracker_recounted_simhits = 1;
669  n_muon_simhits = 1;
670  }
671  n_global_simhits = n_tracker_recounted_simhits + n_muon_simhits;
672 
673  if (UseMuon) {
674  n_muon_selected_simhits = n_muon_simhits;
675  n_global_selected_simhits = n_muon_selected_simhits;
676  }
677  if (UseTracker) {
678  n_tracker_selected_simhits = n_tracker_recounted_simhits;
679  n_global_selected_simhits += n_tracker_selected_simhits;
680  }
681 
683  tracker_quality = static_cast<double>(tracker_nshared);
684  else if (n_tracker_selected_simhits != 0)
685  tracker_quality = static_cast<double>(tracker_nshared) / static_cast<double>(n_tracker_selected_simhits);
686  else
687  tracker_quality = 0;
688 
690  muon_quality = static_cast<double>(muon_nshared);
691  else if (n_muon_selected_simhits != 0)
692  muon_quality = static_cast<double>(muon_nshared) / static_cast<double>(n_muon_selected_simhits);
693  else
694  muon_quality = 0;
695 
696  // global_quality used to order the matching tracks
697  if (n_global_selected_simhits != 0) {
699  global_quality = global_nshared;
700  else
701  global_quality = static_cast<double>(global_nshared) / static_cast<double>(n_global_selected_simhits);
702  } else
703  global_quality = 0;
704 
705  // global purity
706  if (n_selected_hits != 0) {
708  global_purity = global_nshared;
709  else
710  global_purity = static_cast<double>(global_nshared) / static_cast<double>(n_selected_hits);
711  } else
712  global_purity = 0;
713 
714  bool trackerOk = false;
715  if (n_tracker_selected_hits != 0) {
716  if (tracker_quality > tracker_quality_cut)
717  trackerOk = true;
718 
719  tracker_purity = static_cast<double>(tracker_nshared) / static_cast<double>(n_tracker_selected_hits);
721  tracker_purity = static_cast<double>(tracker_nshared);
722 
723  if ((!AbsoluteNumberOfHits_track) && tracker_purity <= PurityCut_track)
724  trackerOk = false;
725 
726  // if a track has just 3 hits in the Tracker we require that all 3
727  // hits are shared
728  if (ThreeHitTracksAreSpecial && n_tracker_selected_hits == 3 && tracker_nshared < 3)
729  trackerOk = false;
730  }
731 
732  bool muonOk = false;
733  if (n_muon_selected_hits != 0) {
734  if (muon_quality > muon_quality_cut)
735  muonOk = true;
736 
737  muon_purity = static_cast<double>(muon_nshared) / static_cast<double>(n_muon_selected_hits);
739  muon_purity = static_cast<double>(muon_nshared);
740 
741  if ((!AbsoluteNumberOfHits_muon) && muon_purity <= PurityCut_muon)
742  muonOk = false;
743  }
744 
745  // (matchOk) has to account for different track types (tracker-only,
746  // standalone muons, global muons)
747  bool matchOk = trackerOk || muonOk;
748 
749  // only for global tracks: match both tracker and muon stub (if acceptOneStubMatchings==false)
750  // depending on the muon selection reject tracks with only one stub (if rejectBadGlobal==true)
751  //
753  ((n_tracker_selected_hits != 0 && n_muon_selected_hits != 0) || rejectBadGlobal))
754  matchOk = trackerOk && muonOk;
755 
756  if (matchOk) {
757  outputCollection[tpindex].push_back(IndexMatch(tindex, global_quality));
758  any_trackingParticle_matched = true;
759 
760  edm::LogVerbatim("MuonAssociatorByHitsHelper")
761  << "*************************************************************"
762  "***********************************************************"
763  << "\n"
764  << "TrackingParticle: q = " << (*trpart).charge() << ", p = " << (*trpart).p()
765  << ", pT = " << (*trpart).pt() << ", eta = " << (*trpart).eta() << ", phi = " << (*trpart).phi() << "\n"
766  << " pdg code = " << (*trpart).pdgId() << ", made of " << (*trpart).numberOfHits()
767  << " PSimHits, recounted " << n_global_simhits << " PSimHits"
768  << " (tracker:" << n_tracker_recounted_simhits << "/muons:" << n_muon_simhits << ")"
769  << ", from " << (*trpart).g4Tracks().size() << " SimTrack:";
770  for (TrackingParticle::g4t_iterator g4T = (*trpart).g4Track_begin(); g4T != (*trpart).g4Track_end(); ++g4T) {
771  edm::LogVerbatim("MuonAssociatorByHitsHelper")
772  << " Id:" << (*g4T).trackId() << "/Evt:(" << (*g4T).eventId().event() << ","
773  << (*g4T).eventId().bunchCrossing() << ")";
774  }
775  edm::LogVerbatim("MuonAssociatorByHitsHelper")
776  << "\t selected " << n_global_selected_simhits << " PSimHits"
777  << " (tracker:" << n_tracker_selected_simhits << "/muons:" << n_muon_selected_simhits << ")"
778  << "\n\t **MATCHED** with quality = " << global_quality << " (tracker: " << tracker_quality
779  << " / muon: " << muon_quality << ")"
780  << "\n\t and purity = " << global_purity << " (tracker: " << tracker_purity
781  << " / muon: " << muon_purity << ")"
782  << "\n\t N shared hits = " << global_nshared << " (tracker: " << tracker_nshared
783  << " / muon: " << muon_nshared << ")"
784  << "\n"
785  << " to: reco::Track " << tindex << ZeroHitMuon << "\n\t"
786  << " made of " << n_selected_hits << " RecHits (tracker:" << n_tracker_valid
787  << "/muons:" << n_muon_selected_hits << ")";
788  } else {
789  // print something only if this TrackingParticle shares some hits with
790  // the current reco::Track
791  if (global_nshared != 0) {
792  if (printRtS)
793  edm::LogVerbatim("MuonAssociatorByHitsHelper")
794  << "*********************************************************"
795  "*********************************************************"
796  "******"
797  << "\n"
798  << "TrackingParticle: q = " << (*trpart).charge() << ", p = " << (*trpart).p()
799  << ", pT = " << (*trpart).pt() << ", eta = " << (*trpart).eta() << ", phi = " << (*trpart).phi()
800  << "\n"
801  << " pdg code = " << (*trpart).pdgId() << ", made of " << (*trpart).numberOfHits()
802  << " PSimHits, recounted " << n_global_simhits << " PSimHits"
803  << " (tracker:" << n_tracker_recounted_simhits << "/muons:" << n_muon_simhits << ")"
804  << ", from " << (*trpart).g4Tracks().size() << " SimTrack:";
805  for (TrackingParticle::g4t_iterator g4T = (*trpart).g4Track_begin(); g4T != (*trpart).g4Track_end();
806  ++g4T) {
807  if (printRtS)
808  edm::LogVerbatim("MuonAssociatorByHitsHelper")
809  << " Id:" << (*g4T).trackId() << "/Evt:(" << (*g4T).eventId().event() << ","
810  << (*g4T).eventId().bunchCrossing() << ")";
811  }
812  if (printRtS)
813  edm::LogVerbatim("MuonAssociatorByHitsHelper")
814  << "\t selected " << n_global_selected_simhits << " PSimHits"
815  << " (tracker:" << n_tracker_selected_simhits << "/muons:" << n_muon_selected_simhits << ")"
816  << "\n\t NOT matched to reco::Track " << tindex << ZeroHitMuon
817  << " with quality = " << global_quality << " (tracker: " << tracker_quality
818  << " / muon: " << muon_quality << ")"
819  << "\n\t and purity = " << global_purity << " (tracker: " << tracker_purity
820  << " / muon: " << muon_purity << ")"
821  << "\n\t N shared hits = " << global_nshared << " (tracker: " << tracker_nshared
822  << " / muon: " << muon_nshared << ")";
823  }
824  }
825  } // loop over TrackingParticle's
826  } // if(n_matching_simhits != 0)
827  } // loop over reco Tracks
828 
829  if (!any_trackingParticle_matched) {
830  edm::LogVerbatim("MuonAssociatorByHitsHelper")
831  << "\n"
832  << "*******************************************************************"
833  "*****************************************************"
834  << "\n NO TrackingParticle associated to ANY input reco::Track ! \n"
835  << "*******************************************************************"
836  "*****************************************************"
837  << "\n";
838  } else {
839  edm::LogVerbatim("MuonAssociatorByHitsHelper")
840  << "*******************************************************************"
841  "*****************************************************"
842  << "\n";
843  }
844 
845  for (IndexAssociation::iterator it = outputCollection.begin(), ed = outputCollection.end(); it != ed; ++it) {
846  std::sort(it->second.begin(), it->second.end());
847  }
848  return outputCollection;
849 }
850 
852  MapOfMatchedIds &muon_matchedIds_valid,
853  MapOfMatchedIds &tracker_matchedIds_INVALID,
854  MapOfMatchedIds &muon_matchedIds_INVALID,
855  int &n_tracker_valid,
856  int &n_dt_valid,
857  int &n_csc_valid,
858  int &n_rpc_valid,
859  int &n_gem_valid,
860  int &n_tracker_matched_valid,
861  int &n_dt_matched_valid,
862  int &n_csc_matched_valid,
863  int &n_rpc_matched_valid,
864  int &n_gem_matched_valid,
865  int &n_tracker_INVALID,
866  int &n_dt_INVALID,
867  int &n_csc_INVALID,
868  int &n_rpc_INVALID,
869  int &n_gem_INVALID,
870  int &n_tracker_matched_INVALID,
871  int &n_dt_matched_INVALID,
872  int &n_csc_matched_INVALID,
873  int &n_rpc_matched_INVALID,
874  int &n_gem_matched_INVALID,
877  const TrackerHitAssociator *trackertruth,
878  const DTHitAssociator &dttruth,
879  const CSCHitAssociator &csctruth,
880  const RPCHitAssociator &rpctruth,
881  const GEMHitAssociator &gemtruth,
882  bool printRtS,
883  const TrackerTopology *tTopo) const
884 
885 {
886  tracker_matchedIds_valid.clear();
887  muon_matchedIds_valid.clear();
888 
889  tracker_matchedIds_INVALID.clear();
890  muon_matchedIds_INVALID.clear();
891 
892  n_tracker_valid = 0;
893  n_dt_valid = 0;
894  n_csc_valid = 0;
895  n_rpc_valid = 0;
896  n_gem_valid = 0;
897 
898  n_tracker_matched_valid = 0;
899  n_dt_matched_valid = 0;
900  n_csc_matched_valid = 0;
901  n_rpc_matched_valid = 0;
902  n_gem_matched_valid = 0;
903 
904  n_tracker_INVALID = 0;
905  n_dt_INVALID = 0;
906  n_csc_INVALID = 0;
907  n_rpc_INVALID = 0;
908  n_gem_INVALID = 0;
909 
910  n_tracker_matched_INVALID = 0;
911  n_dt_matched_INVALID = 0;
912  n_csc_matched_INVALID = 0;
913  n_rpc_matched_INVALID = 0;
914  n_gem_matched_INVALID = 0;
915 
916  std::vector<SimHitIdpr> SimTrackIds;
917 
918  // main loop on TrackingRecHits
919  int iloop = 0;
920  int iH = -1;
921  for (trackingRecHit_iterator it = begin; it != end; it++, iloop++) {
922  stringstream hit_index;
923  hit_index << iloop;
924 
925  const TrackingRecHit *hitp = getHitPtr(it);
926  DetId geoid = hitp->geographicalId();
927 
928  unsigned int detid = geoid.rawId();
929  stringstream detector_id;
930  detector_id << detid;
931 
932  string hitlog = "TrackingRecHit " + hit_index.str();
933  string wireidlog;
934  std::vector<string> DTSimHits;
935 
936  DetId::Detector det = geoid.det();
937  int subdet = geoid.subdetId();
938 
939  bool valid_Hit = hitp->isValid();
940 
941  // Si-Tracker Hits
942  if (det == DetId::Tracker && UseTracker) {
943  stringstream detector_id;
944  detector_id << tTopo->print(detid);
945 
946  if (valid_Hit)
947  hitlog = hitlog + " -Tracker - detID = " + detector_id.str();
948  else
949  hitlog = hitlog + " *** INVALID ***" + " -Tracker - detID = " + detector_id.str();
950 
951  iH++;
952  SimTrackIds = trackertruth->associateHitId(*hitp);
953 
954  if (valid_Hit) {
955  n_tracker_valid++;
956 
957  if (!SimTrackIds.empty()) {
958  n_tracker_matched_valid++;
959  // tracker_matchedIds_valid[iH] = SimTrackIds;
960  tracker_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
961  }
962  } else {
963  n_tracker_INVALID++;
964 
965  if (!SimTrackIds.empty()) {
966  n_tracker_matched_INVALID++;
967  // tracker_matchedIds_INVALID[iH] = SimTrackIds;
968  tracker_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
969  }
970  }
971  }
972  // Muon detector Hits
973  else if (det == DetId::Muon && UseMuon) {
974  // DT Hits
975  if (subdet == MuonSubdetId::DT) {
976  DTWireId dtdetid = DTWireId(detid);
977  stringstream dt_detector_id;
978  dt_detector_id << dtdetid;
979  if (valid_Hit)
980  hitlog = hitlog + " -Muon DT - detID = " + dt_detector_id.str();
981  else
982  hitlog = hitlog + " *** INVALID ***" + " -Muon DT - detID = " + dt_detector_id.str();
983 
984  const DTRecHit1D *dtrechit = dynamic_cast<const DTRecHit1D *>(hitp);
985 
986  // single DT hits
987  if (dtrechit) {
988  iH++;
989  SimTrackIds = dttruth.associateDTHitId(dtrechit);
990 
991  if (valid_Hit) {
992  n_dt_valid++;
993 
994  if (!SimTrackIds.empty()) {
995  n_dt_matched_valid++;
996  // muon_matchedIds_valid[iH] = SimTrackIds;
997  muon_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
998  }
999  } else {
1000  n_dt_INVALID++;
1001 
1002  if (!SimTrackIds.empty()) {
1003  n_dt_matched_INVALID++;
1004  // muon_matchedIds_INVALID[iH] = SimTrackIds;
1005  muon_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
1006  }
1007  }
1008 
1009  if (dumpDT) {
1010  DTWireId wireid = dtrechit->wireId();
1011  stringstream wid;
1012  wid << wireid;
1013  std::vector<PSimHit> dtSimHits = dttruth.associateHit(*hitp);
1014 
1015  stringstream ndthits;
1016  ndthits << dtSimHits.size();
1017  wireidlog = "\t DTWireId :" + wid.str() + ", " + ndthits.str() + " associated PSimHit :";
1018 
1019  for (unsigned int j = 0; j < dtSimHits.size(); j++) {
1020  stringstream index;
1021  index << j;
1022  stringstream simhit;
1023  simhit << dtSimHits[j];
1024  string simhitlog = "\t\t PSimHit " + index.str() + ": " + simhit.str();
1025  DTSimHits.push_back(simhitlog);
1026  }
1027  } // if (dumpDT)
1028  }
1029 
1030  // DT segments
1031  else {
1032  const DTRecSegment4D *dtsegment = dynamic_cast<const DTRecSegment4D *>(hitp);
1033 
1034  if (dtsegment) {
1035  std::vector<const TrackingRecHit *> componentHits, phiHits, zHits;
1036  if (dtsegment->hasPhi()) {
1037  phiHits = dtsegment->phiSegment()->recHits();
1038  componentHits.insert(componentHits.end(), phiHits.begin(), phiHits.end());
1039  }
1040  if (dtsegment->hasZed()) {
1041  zHits = dtsegment->zSegment()->recHits();
1042  componentHits.insert(componentHits.end(), zHits.begin(), zHits.end());
1043  }
1044  if (printRtS)
1045  edm::LogVerbatim("MuonAssociatorByHitsHelper")
1046  << "\n\t this TrackingRecHit is a DTRecSegment4D with " << componentHits.size()
1047  << " hits (phi:" << phiHits.size() << ", z:" << zHits.size() << ")";
1048 
1049  SimTrackIds.clear();
1050  std::vector<SimHitIdpr> i_SimTrackIds;
1051  int i_compHit = 0;
1052  for (std::vector<const TrackingRecHit *>::const_iterator ithit = componentHits.begin();
1053  ithit != componentHits.end();
1054  ++ithit) {
1055  i_compHit++;
1056 
1057  const DTRecHit1D *dtrechit1D = dynamic_cast<const DTRecHit1D *>(*ithit);
1058 
1059  i_SimTrackIds.clear();
1060  if (dtrechit1D) {
1061  iH++;
1062  i_SimTrackIds = dttruth.associateDTHitId(dtrechit1D);
1063 
1064  if (valid_Hit) {
1065  // validity check is on the segment, but hits are counted
1066  // one-by-one
1067  n_dt_valid++;
1068 
1069  if (!i_SimTrackIds.empty()) {
1070  n_dt_matched_valid++;
1071  // muon_matchedIds_valid[iH] = i_SimTrackIds;
1072  muon_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, i_SimTrackIds));
1073  }
1074  } else {
1075  n_dt_INVALID++;
1076 
1077  if (!i_SimTrackIds.empty()) {
1078  n_dt_matched_INVALID++;
1079  // muon_matchedIds_INVALID[iH] = i_SimTrackIds;
1080  muon_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, i_SimTrackIds));
1081  }
1082  }
1083  } else if (printRtS)
1084  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "*** WARNING in "
1085  "MuonAssociatorByHitsHelper::getMatchedIds, null "
1086  "dynamic_cast of a DT TrackingRecHit !";
1087 
1088  unsigned int i_detid = (*ithit)->geographicalId().rawId();
1089  DTWireId i_dtdetid = DTWireId(i_detid);
1090 
1091  stringstream i_dt_detector_id;
1092  i_dt_detector_id << i_dtdetid;
1093 
1094  stringstream i_ss;
1095  i_ss << "\t\t hit " << i_compHit << " -Muon DT - detID = " << i_dt_detector_id.str();
1096 
1097  string i_hitlog = i_ss.str();
1098  i_hitlog = i_hitlog + write_matched_simtracks(i_SimTrackIds);
1099  if (printRtS)
1100  edm::LogVerbatim("MuonAssociatorByHitsHelper") << i_hitlog;
1101 
1102  SimTrackIds.insert(SimTrackIds.end(), i_SimTrackIds.begin(), i_SimTrackIds.end());
1103  }
1104  } // if (dtsegment)
1105 
1106  else if (printRtS)
1107  edm::LogVerbatim("MuonAssociatorByHitsHelper")
1108  << "*** WARNING in MuonAssociatorByHitsHelper::getMatchedIds, "
1109  "DT TrackingRecHit is neither DTRecHit1D nor DTRecSegment4D "
1110  "! ";
1111  }
1112  }
1113 
1114  // CSC Hits
1115  else if (subdet == MuonSubdetId::CSC) {
1116  CSCDetId cscdetid = CSCDetId(detid);
1117  stringstream csc_detector_id;
1118  csc_detector_id << cscdetid;
1119  if (valid_Hit)
1120  hitlog = hitlog + " -Muon CSC- detID = " + csc_detector_id.str();
1121  else
1122  hitlog = hitlog + " *** INVALID ***" + " -Muon CSC- detID = " + csc_detector_id.str();
1123 
1124  const CSCRecHit2D *cscrechit = dynamic_cast<const CSCRecHit2D *>(hitp);
1125 
1126  // single CSC hits
1127  if (cscrechit) {
1128  iH++;
1129  SimTrackIds = csctruth.associateCSCHitId(cscrechit);
1130 
1131  if (valid_Hit) {
1132  n_csc_valid++;
1133 
1134  if (!SimTrackIds.empty()) {
1135  n_csc_matched_valid++;
1136  // muon_matchedIds_valid[iH] = SimTrackIds;
1137  muon_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
1138  }
1139  } else {
1140  n_csc_INVALID++;
1141 
1142  if (!SimTrackIds.empty()) {
1143  n_csc_matched_INVALID++;
1144  // muon_matchedIds_INVALID[iH] = SimTrackIds;
1145  muon_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
1146  }
1147  }
1148  }
1149 
1150  // CSC segments
1151  else {
1152  const CSCSegment *cscsegment = dynamic_cast<const CSCSegment *>(hitp);
1153 
1154  if (cscsegment) {
1155  std::vector<const TrackingRecHit *> componentHits = cscsegment->recHits();
1156  if (printRtS)
1157  edm::LogVerbatim("MuonAssociatorByHitsHelper")
1158  << "\n\t this TrackingRecHit is a CSCSegment with " << componentHits.size() << " hits";
1159 
1160  SimTrackIds.clear();
1161  std::vector<SimHitIdpr> i_SimTrackIds;
1162  int i_compHit = 0;
1163  for (std::vector<const TrackingRecHit *>::const_iterator ithit = componentHits.begin();
1164  ithit != componentHits.end();
1165  ++ithit) {
1166  i_compHit++;
1167 
1168  const CSCRecHit2D *cscrechit2D = dynamic_cast<const CSCRecHit2D *>(*ithit);
1169 
1170  i_SimTrackIds.clear();
1171  if (cscrechit2D) {
1172  iH++;
1173  i_SimTrackIds = csctruth.associateCSCHitId(cscrechit2D);
1174 
1175  if (valid_Hit) {
1176  // validity check is on the segment, but hits are counted
1177  // one-by-one
1178  n_csc_valid++;
1179 
1180  if (!i_SimTrackIds.empty()) {
1181  n_csc_matched_valid++;
1182  // muon_matchedIds_valid[iH] = i_SimTrackIds;
1183  muon_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, i_SimTrackIds));
1184  }
1185  } else {
1186  n_csc_INVALID++;
1187 
1188  if (!i_SimTrackIds.empty()) {
1189  n_csc_matched_INVALID++;
1190  // muon_matchedIds_INVALID[iH] = i_SimTrackIds;
1191  muon_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, i_SimTrackIds));
1192  }
1193  }
1194  } else if (printRtS)
1195  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "*** WARNING in "
1196  "MuonAssociatorByHitsHelper::getMatchedIds, null "
1197  "dynamic_cast of a CSC TrackingRecHit !";
1198 
1199  unsigned int i_detid = (*ithit)->geographicalId().rawId();
1200  CSCDetId i_cscdetid = CSCDetId(i_detid);
1201 
1202  stringstream i_csc_detector_id;
1203  i_csc_detector_id << i_cscdetid;
1204 
1205  stringstream i_ss;
1206  i_ss << "\t\t hit " << i_compHit << " -Muon CSC- detID = " << i_csc_detector_id.str();
1207 
1208  string i_hitlog = i_ss.str();
1209  i_hitlog = i_hitlog + write_matched_simtracks(i_SimTrackIds);
1210  if (printRtS)
1211  edm::LogVerbatim("MuonAssociatorByHitsHelper") << i_hitlog;
1212 
1213  SimTrackIds.insert(SimTrackIds.end(), i_SimTrackIds.begin(), i_SimTrackIds.end());
1214  }
1215  } // if (cscsegment)
1216 
1217  else if (printRtS)
1218  edm::LogVerbatim("MuonAssociatorByHitsHelper")
1219  << "*** WARNING in MuonAssociatorByHitsHelper::getMatchedIds, "
1220  "CSC TrackingRecHit is neither CSCRecHit2D nor CSCSegment "
1221  "! ";
1222  }
1223  }
1224 
1225  // RPC Hits
1226  else if (subdet == MuonSubdetId::RPC) {
1227  RPCDetId rpcdetid = RPCDetId(detid);
1228  stringstream rpc_detector_id;
1229  rpc_detector_id << rpcdetid;
1230  if (valid_Hit)
1231  hitlog = hitlog + " -Muon RPC- detID = " + rpc_detector_id.str();
1232  else
1233  hitlog = hitlog + " *** INVALID ***" + " -Muon RPC- detID = " + rpc_detector_id.str();
1234 
1235  iH++;
1236  SimTrackIds = rpctruth.associateRecHit(*hitp);
1237 
1238  if (valid_Hit) {
1239  n_rpc_valid++;
1240 
1241  if (!SimTrackIds.empty()) {
1242  n_rpc_matched_valid++;
1243  // muon_matchedIds_valid[iH] = SimTrackIds;
1244  muon_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
1245  }
1246  } else {
1247  n_rpc_INVALID++;
1248 
1249  if (!SimTrackIds.empty()) {
1250  n_rpc_matched_INVALID++;
1251  // muon_matchedIds_INVALID[iH] = SimTrackIds;
1252  muon_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
1253  }
1254  }
1255  }
1256 
1257  // GEM Hits
1258  else if (subdet == MuonSubdetId::GEM) {
1259  GEMDetId gemdetid = GEMDetId(detid);
1260  stringstream gem_detector_id;
1261  gem_detector_id << gemdetid;
1262  if (valid_Hit)
1263  hitlog = hitlog + " -Muon GEM- detID = " + gem_detector_id.str();
1264  else
1265  hitlog = hitlog + " *** INVALID ***" + " -Muon GEM- detID = " + gem_detector_id.str();
1266 
1267  const GEMRecHit *gemrechit = dynamic_cast<const GEMRecHit *>(hitp);
1268  if (gemrechit) {
1269  iH++;
1270  SimTrackIds = gemtruth.associateRecHit(gemrechit);
1271 
1272  if (valid_Hit) {
1273  n_gem_valid++;
1274 
1275  if (!SimTrackIds.empty()) {
1276  n_gem_matched_valid++;
1277  // muon_matchedIds_valid[iH] = SimTrackIds;
1278  muon_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
1279  }
1280  } else {
1281  n_gem_INVALID++;
1282 
1283  if (!SimTrackIds.empty()) {
1284  n_gem_matched_INVALID++;
1285  // muon_matchedIds_INVALID[iH] = SimTrackIds;
1286  muon_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
1287  }
1288  }
1289  } else {
1290  const GEMSegment *gemsegment = dynamic_cast<const GEMSegment *>(hitp);
1291  if (gemsegment) {
1292  std::vector<const TrackingRecHit *> componentHits = gemsegment->recHits();
1293  if (printRtS)
1294  edm::LogVerbatim("MuonAssociatorByHitsHelper")
1295  << "\n\t this TrackingRecHit is a GEMSegment with " << componentHits.size() << " hits";
1296 
1297  SimTrackIds.clear();
1298  std::vector<SimHitIdpr> i_SimTrackIds;
1299  int i_compHit = 0;
1300 
1301  for (auto const &ithit : componentHits) {
1302  i_compHit++;
1303 
1304  const GEMRecHit *gemrechitseg = dynamic_cast<const GEMRecHit *>(ithit);
1305 
1306  i_SimTrackIds.clear();
1307  if (gemrechitseg) {
1308  iH++;
1309  i_SimTrackIds = gemtruth.associateRecHit(gemrechitseg);
1310 
1311  if (valid_Hit) {
1312  // validity check is on the segment, but hits are counted
1313  // one-by-one
1314  n_gem_valid++;
1315 
1316  if (!i_SimTrackIds.empty()) {
1317  n_gem_matched_valid++;
1318  // muon_matchedIds_valid[iH] = i_SimTrackIds;
1319  muon_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, i_SimTrackIds));
1320  }
1321  } else {
1322  n_gem_INVALID++;
1323 
1324  if (!i_SimTrackIds.empty()) {
1325  n_gem_matched_INVALID++;
1326  // muon_matchedIds_INVALID[iH] = i_SimTrackIds;
1327  muon_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, i_SimTrackIds));
1328  }
1329  }
1330  } else if (printRtS)
1331  edm::LogVerbatim("MuonAssociatorByHitsHelper") << "*** WARNING in "
1332  "MuonAssociatorByHitsHelper::getMatchedIds, null "
1333  "dynamic_cast of a GEM TrackingRecHit !";
1334 
1335  if (printRtS) {
1336  unsigned int i_detid = ithit->geographicalId().rawId();
1337  GEMDetId i_gemdetid = GEMDetId(i_detid);
1338 
1339  string i_hitlog = std::to_string(i_gemdetid);
1340  i_hitlog = i_hitlog + write_matched_simtracks(i_SimTrackIds);
1341  edm::LogVerbatim("MuonAssociatorByHitsHelper") << i_hitlog;
1342  }
1343 
1344  SimTrackIds.insert(SimTrackIds.end(), i_SimTrackIds.begin(), i_SimTrackIds.end());
1345  }
1346  } // if (gemsegment)
1347 
1348  } // if (gemrechit
1349  } else if (printRtS)
1350  edm::LogVerbatim("MuonAssociatorByHitsHelper")
1351  << "TrackingRecHit " << iloop << " *** WARNING *** Unexpected Hit from Detector = " << det;
1352  } // end if (det == DetId::Muon && UseMuon)
1353  else
1354  continue;
1355 
1356  hitlog = hitlog + write_matched_simtracks(SimTrackIds);
1357 
1358  if (printRtS)
1359  edm::LogVerbatim("MuonAssociatorByHitsHelper") << hitlog;
1360  if (printRtS && dumpDT && det == DetId::Muon && subdet == MuonSubdetId::DT) {
1361  edm::LogVerbatim("MuonAssociatorByHitsHelper") << wireidlog;
1362  for (unsigned int j = 0; j < DTSimHits.size(); j++) {
1363  edm::LogVerbatim("MuonAssociatorByHitsHelper") << DTSimHits[j];
1364  }
1365  }
1366 
1367  } // trackingRecHit loop
1368 }
1369 
1371  TrackingParticleCollection::const_iterator trpart) const {
1372  int nshared = 0;
1373  const std::vector<SimTrack> &g4Tracks = trpart->g4Tracks();
1374 
1375  // map is indexed over the rechits of the reco::Track (no double-countings
1376  // allowed)
1377  for (MapOfMatchedIds::const_iterator iRecH = matchedIds.begin(); iRecH != matchedIds.end(); ++iRecH) {
1378  // vector of associated simhits associated to the current rechit
1379  std::vector<SimHitIdpr> const &SimTrackIds = (*iRecH)->second;
1380 
1381  bool found = false;
1382 
1383  for (const auto &iSimH : SimTrackIds) {
1384  uint32_t simtrackId = iSimH.first;
1385  EncodedEventId evtId = iSimH.second;
1386 
1387  // look for shared hits with the given TrackingParticle (looping over
1388  // component SimTracks)
1389  for (const auto &simtrack : g4Tracks) {
1390  if (simtrack.trackId() == simtrackId && simtrack.eventId() == evtId) {
1391  found = true;
1392  break;
1393  }
1394  }
1395 
1396  if (found) {
1397  nshared++;
1398  break;
1399  }
1400  }
1401  }
1402 
1403  return nshared;
1404 }
1405 
1406 std::string MuonAssociatorByHitsHelper::write_matched_simtracks(const std::vector<SimHitIdpr> &SimTrackIds) const {
1407  if (SimTrackIds.empty())
1408  return " *** UNMATCHED ***";
1409 
1410  string hitlog(" matched to SimTrack");
1411 
1412  for (size_t j = 0; j < SimTrackIds.size(); j++) {
1413  char buf[64];
1414  snprintf(buf,
1415  64,
1416  " Id:%i/Evt:(%i,%i) ",
1417  SimTrackIds[j].first,
1418  SimTrackIds[j].second.event(),
1419  SimTrackIds[j].second.bunchCrossing());
1420  hitlog += buf;
1421  }
1422  return hitlog;
1423 }
Log< level::Info, true > LogVerbatim
std::vector< const TrackingRecHit * > recHits() const override
Access to component RecHits (if any)
Definition: CSCSegment.cc:32
static constexpr int GEM
Definition: MuonSubdetId.h:14
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:57
tuple resources
Definition: mps_fire.py:263
std::string to_string(const V &value)
Definition: OMSAccess.h:71
std::vector< SimHitIdpr > associateCSCHitId(const CSCRecHit2D *) const
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
#define LogTrace(id)
U second(std::pair< T, U > const &p)
std::vector< SimHitIdpr > associateRecHit(const TrackingRecHit &hit) const
std::vector< std::unique_ptr< uint_SimHitIdpr_pair > > MapOfMatchedIds
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:48
std::vector< const TrackingRecHit * > recHits() const override
Access to component RecHits (if any)
bool hasPhi() const
Does it have the Phi projection?
std::vector< const TrackingRecHit * > recHits() const override
Access to component RecHits (if any)
Definition: GEMSegment.cc:82
std::vector< SimTrack >::const_iterator g4t_iterator
bool hasZed() const
Does it have the Z projection?
Definition: DetId.h:17
const DTSLRecSegment2D * zSegment() const
The Z segment: 0 if not zed projection available.
std::vector< SimHitIdpr > associateRecHit(const GEMRecHit *gemrechit) const
Detector
Definition: DetId.h:24
static constexpr int RPC
Definition: MuonSubdetId.h:13
bool isValid() const
IndexAssociation associateRecoToSimIndices(const TrackHitsCollection &, const edm::RefVector< TrackingParticleCollection > &, Resources const &) const
std::vector< SimHitIdpr > associateHitId(const TrackingRecHit &thit) const
const TrackingRecHit * getHitPtr(edm::OwnVector< TrackingRecHit >::const_iterator iter) const
string end
Definition: dataset.py:937
size_type size() const
Size of the RefVector.
Definition: RefVector.h:102
std::vector< TrackingParticle > TrackingParticleCollection
DetId geographicalId() const
static constexpr int DT
Definition: MuonSubdetId.h:11
static constexpr int CSC
Definition: MuonSubdetId.h:12
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:46
DTWireId wireId() const
Return the wireId.
Definition: DTRecHit1D.h:76