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