CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonAssociatorByHitsHelper.cc
Go to the documentation of this file.
1 #include "SimMuon/MCTruth/interface/MuonAssociatorByHitsHelper.h"
2 #include "FWCore/MessageLogger/interface/MessageLogger.h"
3 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
4 #include "DataFormats/DetId/interface/DetId.h"
5 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
6 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
7 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
8 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
9 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
10 #include "DataFormats/DTRecHit/interface/DTRecSegment4D.h"
11 #include "DataFormats/CSCRecHit/interface/CSCSegment.h"
12 #include <sstream>
13 
14 using namespace reco;
15 using namespace std;
16 
17 MuonAssociatorByHitsHelper::MuonAssociatorByHitsHelper (const edm::ParameterSet& conf) :
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 
58 MuonAssociatorByHitsHelper::IndexAssociation
59 MuonAssociatorByHitsHelper::associateRecoToSimIndices(const TrackHitsCollection & tC,
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 
92  TrackingParticleCollection tPC;
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::LogWarning("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) {
263  if (AbsoluteNumberOfHits_muon && AbsoluteNumberOfHits_track)
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.size())
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 
340 MuonAssociatorByHitsHelper::IndexAssociation
341 MuonAssociatorByHitsHelper::associateSimToRecoIndices( const TrackHitsCollection & tC,
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 
379  TrackingParticleCollection tPC;
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::LogWarning("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  }
620  n_global_simhits = n_tracker_recounted_simhits + n_muon_simhits;
621 
622  if (UseMuon) {
623  n_muon_selected_simhits = n_muon_simhits;
624  n_global_selected_simhits = n_muon_selected_simhits;
625  }
626  if (UseTracker) {
627  n_tracker_selected_simhits = n_tracker_recounted_simhits;
628  n_global_selected_simhits += n_tracker_selected_simhits;
629  }
630 
631  if (AbsoluteNumberOfHits_track) tracker_quality = static_cast<double>(tracker_nshared);
632  else if (n_tracker_selected_simhits!=0)
633  tracker_quality = static_cast<double>(tracker_nshared)/static_cast<double>(n_tracker_selected_simhits);
634  else tracker_quality = 0;
635 
636  if (AbsoluteNumberOfHits_muon) muon_quality = static_cast<double>(muon_nshared);
637  else if (n_muon_selected_simhits!=0)
638  muon_quality = static_cast<double>(muon_nshared)/static_cast<double>(n_muon_selected_simhits);
639  else muon_quality = 0;
640 
641  // global_quality used to order the matching tracks
642  if (n_global_selected_simhits != 0) {
643  if (AbsoluteNumberOfHits_muon && AbsoluteNumberOfHits_track)
644  global_quality = global_nshared;
645  else
646  global_quality = static_cast<double>(global_nshared)/static_cast<double>(n_global_selected_simhits);
647  }
648  else global_quality = 0;
649 
650  // global purity
651  if (n_selected_hits != 0) {
652  if (AbsoluteNumberOfHits_muon && AbsoluteNumberOfHits_track)
653  global_purity = global_nshared;
654  else
655  global_purity = static_cast<double>(global_nshared)/static_cast<double>(n_selected_hits);
656  }
657  else global_purity = 0;
658 
659  bool trackerOk = false;
660  if (n_tracker_selected_hits != 0) {
661  if (tracker_quality > tracker_quality_cut) trackerOk = true;
662 
663  tracker_purity = static_cast<double>(tracker_nshared)/static_cast<double>(n_tracker_selected_hits);
664  if (AbsoluteNumberOfHits_track) tracker_purity = static_cast<double>(tracker_nshared);
665 
666  if ((!AbsoluteNumberOfHits_track) && tracker_purity <= PurityCut_track) trackerOk = false;
667 
668  //if a track has just 3 hits in the Tracker we require that all 3 hits are shared
669  if (ThreeHitTracksAreSpecial && n_tracker_selected_hits==3 && tracker_nshared<3) trackerOk = false;
670  }
671 
672  bool muonOk = false;
673  if (n_muon_selected_hits != 0) {
674  if (muon_quality > muon_quality_cut) muonOk = true;
675 
676  muon_purity = static_cast<double>(muon_nshared)/static_cast<double>(n_muon_selected_hits);
677  if (AbsoluteNumberOfHits_muon) muon_purity = static_cast<double>(muon_nshared);
678 
679  if ((!AbsoluteNumberOfHits_muon) && muon_purity <= PurityCut_muon) muonOk = false;
680  }
681 
682  // (matchOk) has to account for different track types (tracker-only, standalone muons, global muons)
683  bool matchOk = trackerOk || muonOk;
684 
685  // only for global muons: match both tracker and muon stub unless (acceptOneStubMatchings==true)
686  if (!acceptOneStubMatchings && n_tracker_selected_hits!=0 && n_muon_selected_hits!=0)
687  matchOk = trackerOk && muonOk;
688 
689  if (matchOk) {
690 
691  outputCollection[tpindex].push_back(IndexMatch(tindex,global_quality));
692  any_trackingParticle_matched = true;
693 
694  edm::LogVerbatim("MuonAssociatorByHitsHelper")
695  <<"************************************************************************************************************************"
696  <<"\n"<< "TrackingParticle " << tpindex <<", q = "<<(*trpart).charge()<<", p = "<<(*trpart).p()
697  <<", pT = "<<(*trpart).pt()<<", eta = "<<(*trpart).eta()<<", phi = "<<(*trpart).phi()
698  <<"\n"<<" pdg code = "<<(*trpart).pdgId()
699  <<", made of "<<(*trpart).numberOfHits()<<" PSimHits, recounted "<<n_global_simhits<<" PSimHits"
700  <<" (tracker:"<<n_tracker_recounted_simhits<<"/muons:"<<n_muon_simhits<<")"
701  <<", from "<<(*trpart).g4Tracks().size()<<" SimTrack:";
702  for(TrackingParticle::g4t_iterator g4T=(*trpart).g4Track_begin();
703  g4T!=(*trpart).g4Track_end();
704  ++g4T) {
705  edm::LogVerbatim("MuonAssociatorByHitsHelper")
706  <<" Id:"<<(*g4T).trackId()<<"/Evt:("<<(*g4T).eventId().event()<<","<<(*g4T).eventId().bunchCrossing()<<")";
707  }
708  edm::LogVerbatim("MuonAssociatorByHitsHelper")
709  <<"\t selected "<<n_global_selected_simhits<<" PSimHits"
710  <<" (tracker:"<<n_tracker_selected_simhits<<"/muons:"<<n_muon_selected_simhits<<")"
711  << "\n\t **MATCHED** with quality = "<<global_quality<<" (tracker: "<<tracker_quality<<" / muon: "<<muon_quality<<")"
712  << "\n\t and purity = "<<global_purity<<" (tracker: "<<tracker_purity<<" / muon: "<<muon_purity<<")"
713  << "\n\t N shared hits = "<<global_nshared<<" (tracker: "<<tracker_nshared<<" / muon: "<<muon_nshared<<")"
714  <<"\n" <<" to: reco::Track "<<tindex<<ZeroHitMuon
715  <<"\n\t"<< " made of "<<n_selected_hits<<" RecHits (tracker:"<<n_tracker_valid<<"/muons:"<<n_muon_selected_hits<<")";
716  }
717  else {
718  // print something only if this TrackingParticle shares some hits with the current reco::Track
719  if (global_nshared != 0) {
720  if (printRtS) edm::LogVerbatim("MuonAssociatorByHitsHelper")
721  <<"************************************************************************************************************************"
722  <<"\n"<<"TrackingParticle " << tpindex <<", q = "<<(*trpart).charge()<<", p = "<<(*trpart).p()
723  <<", pT = "<<(*trpart).pt()<<", eta = "<<(*trpart).eta()<<", phi = "<<(*trpart).phi()
724  <<"\n"<<" pdg code = "<<(*trpart).pdgId()
725  <<", made of "<<(*trpart).numberOfHits()<<" PSimHits, recounted "<<n_global_simhits<<" PSimHits"
726  <<" (tracker:"<<n_tracker_recounted_simhits<<"/muons:"<<n_muon_simhits<<")"
727  <<", from "<<(*trpart).g4Tracks().size()<<" SimTrack:";
728  for(TrackingParticle::g4t_iterator g4T=(*trpart).g4Track_begin();
729  g4T!=(*trpart).g4Track_end();
730  ++g4T) {
731  if (printRtS) edm::LogVerbatim("MuonAssociatorByHitsHelper")
732  <<" Id:"<<(*g4T).trackId()<<"/Evt:("<<(*g4T).eventId().event()<<","<<(*g4T).eventId().bunchCrossing()<<")";
733  }
734  if (printRtS) edm::LogVerbatim("MuonAssociatorByHitsHelper")
735  <<"\t selected "<<n_global_selected_simhits<<" PSimHits"
736  <<" (tracker:"<<n_tracker_selected_simhits<<"/muons:"<<n_muon_selected_simhits<<")"
737  <<"\n\t NOT matched to reco::Track "<<tindex<<ZeroHitMuon
738  <<" with quality = "<<global_quality<<" (tracker: "<<tracker_quality<<" / muon: "<<muon_quality<<")"
739  <<"\n\t and purity = "<<global_purity<<" (tracker: "<<tracker_purity<<" / muon: "<<muon_purity<<")"
740  <<"\n\t N shared hits = "<<global_nshared<<" (tracker: "<<tracker_nshared<<" / muon: "<<muon_nshared<<")";
741  }
742  }
743  } // loop over TrackingParticle's
744  } // if(n_matching_simhits != 0)
745  } // loop over reco Tracks
746 
747  if (!any_trackingParticle_matched) {
748  edm::LogVerbatim("MuonAssociatorByHitsHelper")
749  <<"\n"
750  <<"************************************************************************************************************************"
751  << "\n NO TrackingParticle associated to ANY input reco::Track ! \n"
752  <<"************************************************************************************************************************"<<"\n";
753  } else {
754  edm::LogVerbatim("MuonAssociatorByHitsHelper")
755  <<"************************************************************************************************************************"<<"\n";
756  }
757 
758  for (IndexAssociation::iterator it = outputCollection.begin(), ed = outputCollection.end(); it != ed; ++it) {
759  std::sort(it->second.begin(), it->second.end());
760  }
761  return outputCollection;
762 }
763 
764 
765 void MuonAssociatorByHitsHelper::getMatchedIds
766 (MapOfMatchedIds & tracker_matchedIds_valid, MapOfMatchedIds & muon_matchedIds_valid,
767  MapOfMatchedIds & tracker_matchedIds_INVALID, MapOfMatchedIds & muon_matchedIds_INVALID,
768  int& n_tracker_valid, int& n_dt_valid, int& n_csc_valid, int& n_rpc_valid, int& n_gem_valid,
769  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,
770  int& n_tracker_INVALID, int& n_dt_INVALID, int& n_csc_INVALID, int& n_rpc_INVALID, int& n_gem_INVALID,
771  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,
772  trackingRecHit_iterator begin, trackingRecHit_iterator end,
773  const TrackerHitAssociator* trackertruth,
774  const DTHitAssociator& dttruth, const CSCHitAssociator& csctruth, const RPCHitAssociator& rpctruth, const GEMHitAssociator& gemtruth, bool printRtS,
775  const TrackerTopology *tTopo) const
776 
777 {
778  tracker_matchedIds_valid.clear();
779  muon_matchedIds_valid.clear();
780 
781  tracker_matchedIds_INVALID.clear();
782  muon_matchedIds_INVALID.clear();
783 
784  n_tracker_valid = 0;
785  n_dt_valid = 0;
786  n_csc_valid = 0;
787  n_rpc_valid = 0;
788  n_gem_valid = 0;
789 
790  n_tracker_matched_valid = 0;
791  n_dt_matched_valid = 0;
792  n_csc_matched_valid = 0;
793  n_rpc_matched_valid = 0;
794  n_gem_matched_valid = 0;
795 
796  n_tracker_INVALID = 0;
797  n_dt_INVALID = 0;
798  n_csc_INVALID = 0;
799  n_rpc_INVALID = 0;
800  n_gem_INVALID = 0;
801 
802  n_tracker_matched_INVALID = 0;
803  n_dt_matched_INVALID = 0;
804  n_csc_matched_INVALID = 0;
805  n_rpc_matched_INVALID = 0;
806  n_gem_matched_INVALID = 0;
807 
808  std::vector<SimHitIdpr> SimTrackIds;
809 
810  // main loop on TrackingRecHits
811  int iloop = 0;
812  int iH = -1;
813  for (trackingRecHit_iterator it = begin; it != end; it++, iloop++) {
814  stringstream hit_index;
815  hit_index<<iloop;
816 
817  const TrackingRecHit * hitp = getHitPtr(it);
818  DetId geoid = hitp->geographicalId();
819 
820  unsigned int detid = geoid.rawId();
821  stringstream detector_id;
822  detector_id<<detid;
823 
824  string hitlog = "TrackingRecHit "+hit_index.str();
825  string wireidlog;
826  std::vector<string> DTSimHits;
827 
828  DetId::Detector det = geoid.det();
829  int subdet = geoid.subdetId();
830 
831  bool valid_Hit = hitp->isValid();
832 
833  // Si-Tracker Hits
834  if (det == DetId::Tracker && UseTracker) {
835  stringstream detector_id;
836  detector_id<< tTopo->print(detid);
837 
838  if (valid_Hit) hitlog = hitlog+" -Tracker - detID = "+detector_id.str();
839  else hitlog = hitlog+" *** INVALID ***"+" -Tracker - detID = "+detector_id.str();
840 
841  iH++;
842  SimTrackIds = trackertruth->associateHitId(*hitp);
843 
844  if (valid_Hit) {
845  n_tracker_valid++;
846 
847  if(!SimTrackIds.empty()) {
848  n_tracker_matched_valid++;
849  //tracker_matchedIds_valid[iH] = SimTrackIds;
850  tracker_matchedIds_valid.push_back( new uint_SimHitIdpr_pair(iH, SimTrackIds));
851  }
852  } else {
853  n_tracker_INVALID++;
854 
855  if(!SimTrackIds.empty()) {
856  n_tracker_matched_INVALID++;
857  //tracker_matchedIds_INVALID[iH] = SimTrackIds;
858  tracker_matchedIds_INVALID.push_back( new uint_SimHitIdpr_pair(iH, SimTrackIds));
859  }
860  }
861  }
862  // Muon detector Hits
863  else if (det == DetId::Muon && UseMuon) {
864 
865  // DT Hits
866  if (subdet == MuonSubdetId::DT) {
867  DTWireId dtdetid = DTWireId(detid);
868  stringstream dt_detector_id;
869  dt_detector_id << dtdetid;
870  if (valid_Hit) hitlog = hitlog+" -Muon DT - detID = "+dt_detector_id.str();
871  else hitlog = hitlog+" *** INVALID ***"+" -Muon DT - detID = "+dt_detector_id.str();
872 
873  const DTRecHit1D * dtrechit = dynamic_cast<const DTRecHit1D *>(hitp);
874 
875  // single DT hits
876  if (dtrechit) {
877  iH++;
878  SimTrackIds = dttruth.associateDTHitId(dtrechit);
879 
880  if (valid_Hit) {
881  n_dt_valid++;
882 
883  if (!SimTrackIds.empty()) {
884  n_dt_matched_valid++;
885  //muon_matchedIds_valid[iH] = SimTrackIds;
886  muon_matchedIds_valid.push_back (new uint_SimHitIdpr_pair(iH,SimTrackIds));
887  }
888  } else {
889  n_dt_INVALID++;
890 
891  if (!SimTrackIds.empty()) {
892  n_dt_matched_INVALID++;
893  //muon_matchedIds_INVALID[iH] = SimTrackIds;
894  muon_matchedIds_INVALID.push_back (new uint_SimHitIdpr_pair(iH,SimTrackIds));
895  }
896  }
897 
898  if (dumpDT) {
899  DTWireId wireid = dtrechit->wireId();
900  stringstream wid;
901  wid<<wireid;
902  std::vector<PSimHit> dtSimHits = dttruth.associateHit(*hitp);
903 
904  stringstream ndthits;
905  ndthits<<dtSimHits.size();
906  wireidlog = "\t DTWireId :"+wid.str()+", "+ndthits.str()+" associated PSimHit :";
907 
908  for (unsigned int j=0; j<dtSimHits.size(); j++) {
909  stringstream index;
910  index<<j;
911  stringstream simhit;
912  simhit<<dtSimHits[j];
913  string simhitlog = "\t\t PSimHit "+index.str()+": "+simhit.str();
914  DTSimHits.push_back(simhitlog);
915  }
916  } // if (dumpDT)
917  }
918 
919  // DT segments
920  else {
921  const DTRecSegment4D * dtsegment = dynamic_cast<const DTRecSegment4D *>(hitp);
922 
923  if (dtsegment) {
924 
925  std::vector<const TrackingRecHit *> componentHits, phiHits, zHits;
926  if (dtsegment->hasPhi()) {
927  phiHits = dtsegment->phiSegment()->recHits();
928  componentHits.insert(componentHits.end(),phiHits.begin(),phiHits.end());
929  }
930  if (dtsegment->hasZed()) {
931  zHits = dtsegment->zSegment()->recHits();
932  componentHits.insert(componentHits.end(),zHits.begin(),zHits.end());
933  }
934  if (printRtS) edm::LogVerbatim("MuonAssociatorByHitsHelper")
935  <<"\n\t this TrackingRecHit is a DTRecSegment4D with "
936  <<componentHits.size()<<" hits (phi:"<<phiHits.size()<<", z:"<<zHits.size()<<")";
937 
938  std::vector<SimHitIdpr> i_SimTrackIds;
939  int i_compHit = 0;
940  for (std::vector<const TrackingRecHit *>::const_iterator ithit =componentHits.begin();
941  ithit != componentHits.end(); ++ithit) {
942  i_compHit++;
943 
944  const DTRecHit1D * dtrechit1D = dynamic_cast<const DTRecHit1D *>(*ithit);
945 
946  i_SimTrackIds.clear();
947  if (dtrechit1D) {
948  iH++;
949  i_SimTrackIds = dttruth.associateDTHitId(dtrechit1D);
950 
951  if (valid_Hit) {
952  // validity check is on the segment, but hits are counted one-by-one
953  n_dt_valid++;
954 
955  if (!i_SimTrackIds.empty()) {
956  n_dt_matched_valid++;
957  //muon_matchedIds_valid[iH] = i_SimTrackIds;
958  muon_matchedIds_valid.push_back (new uint_SimHitIdpr_pair(iH,i_SimTrackIds));
959  }
960  } else {
961  n_dt_INVALID++;
962 
963  if (!i_SimTrackIds.empty()) {
964  n_dt_matched_INVALID++;
965  //muon_matchedIds_INVALID[iH] = i_SimTrackIds;
966  muon_matchedIds_INVALID.push_back (new uint_SimHitIdpr_pair(iH,i_SimTrackIds));
967 
968  }
969  }
970  } else if (printRtS) edm::LogWarning("MuonAssociatorByHitsHelper")
971  <<"*** WARNING in MuonAssociatorByHitsHelper::getMatchedIds, null dynamic_cast of a DT TrackingRecHit !";
972 
973  unsigned int i_detid = (*ithit)->geographicalId().rawId();
974  DTWireId i_dtdetid = DTWireId(i_detid);
975 
976  stringstream i_dt_detector_id;
977  i_dt_detector_id << i_dtdetid;
978 
979  stringstream i_ss;
980  i_ss<<"\t\t hit "<<i_compHit<<" -Muon DT - detID = "<<i_dt_detector_id.str();
981 
982  string i_hitlog = i_ss.str();
983  i_hitlog = i_hitlog + write_matched_simtracks(i_SimTrackIds);
984  if (printRtS) edm::LogVerbatim("MuonAssociatorByHitsHelper") << i_hitlog;
985 
986  SimTrackIds.insert(SimTrackIds.end(),i_SimTrackIds.begin(),i_SimTrackIds.end());
987  }
988  } // if (dtsegment)
989 
990  else if (printRtS) edm::LogWarning("MuonAssociatorByHitsHelper")
991  <<"*** WARNING in MuonAssociatorByHitsHelper::getMatchedIds, DT TrackingRecHit is neither DTRecHit1D nor DTRecSegment4D ! ";
992  }
993  }
994 
995  // CSC Hits
996  else if (subdet == MuonSubdetId::CSC) {
997  CSCDetId cscdetid = CSCDetId(detid);
998  stringstream csc_detector_id;
999  csc_detector_id << cscdetid;
1000  if (valid_Hit) hitlog = hitlog+" -Muon CSC- detID = "+csc_detector_id.str();
1001  else hitlog = hitlog+" *** INVALID ***"+" -Muon CSC- detID = "+csc_detector_id.str();
1002 
1003  const CSCRecHit2D * cscrechit = dynamic_cast<const CSCRecHit2D *>(hitp);
1004 
1005  // single CSC hits
1006  if (cscrechit) {
1007  iH++;
1008  SimTrackIds = csctruth.associateCSCHitId(cscrechit);
1009 
1010  if (valid_Hit) {
1011  n_csc_valid++;
1012 
1013  if (!SimTrackIds.empty()) {
1014  n_csc_matched_valid++;
1015  //muon_matchedIds_valid[iH] = SimTrackIds;
1016  muon_matchedIds_valid.push_back (new uint_SimHitIdpr_pair(iH,SimTrackIds));
1017  }
1018  } else {
1019  n_csc_INVALID++;
1020 
1021  if (!SimTrackIds.empty()) {
1022  n_csc_matched_INVALID++;
1023  //muon_matchedIds_INVALID[iH] = SimTrackIds;
1024  muon_matchedIds_INVALID.push_back (new uint_SimHitIdpr_pair(iH,SimTrackIds));
1025  }
1026  }
1027  }
1028 
1029  // CSC segments
1030  else {
1031  const CSCSegment * cscsegment = dynamic_cast<const CSCSegment *>(hitp);
1032 
1033  if (cscsegment) {
1034 
1035  std::vector<const TrackingRecHit *> componentHits = cscsegment->recHits();
1036  if (printRtS) edm::LogVerbatim("MuonAssociatorByHitsHelper")
1037  <<"\n\t this TrackingRecHit is a CSCSegment with "<<componentHits.size()<<" hits";
1038 
1039  std::vector<SimHitIdpr> i_SimTrackIds;
1040  int i_compHit = 0;
1041  for (std::vector<const TrackingRecHit *>::const_iterator ithit =componentHits.begin();
1042  ithit != componentHits.end(); ++ithit) {
1043  i_compHit++;
1044 
1045  const CSCRecHit2D * cscrechit2D = dynamic_cast<const CSCRecHit2D *>(*ithit);
1046 
1047  i_SimTrackIds.clear();
1048  if (cscrechit2D) {
1049  iH++;
1050  i_SimTrackIds = csctruth.associateCSCHitId(cscrechit2D);
1051 
1052  if (valid_Hit) {
1053  // validity check is on the segment, but hits are counted one-by-one
1054  n_csc_valid++;
1055 
1056  if (!i_SimTrackIds.empty()) {
1057  n_csc_matched_valid++;
1058  //muon_matchedIds_valid[iH] = i_SimTrackIds;
1059  muon_matchedIds_valid.push_back (new uint_SimHitIdpr_pair(iH,i_SimTrackIds));
1060  }
1061  } else {
1062  n_csc_INVALID++;
1063 
1064  if (!i_SimTrackIds.empty()) {
1065  n_csc_matched_INVALID++;
1066  //muon_matchedIds_INVALID[iH] = i_SimTrackIds;
1067  muon_matchedIds_INVALID.push_back (new uint_SimHitIdpr_pair(iH,i_SimTrackIds));
1068  }
1069  }
1070  } else if (printRtS) edm::LogWarning("MuonAssociatorByHitsHelper")
1071  <<"*** WARNING in MuonAssociatorByHitsHelper::getMatchedIds, null dynamic_cast of a CSC TrackingRecHit !";
1072 
1073  unsigned int i_detid = (*ithit)->geographicalId().rawId();
1074  CSCDetId i_cscdetid = CSCDetId(i_detid);
1075 
1076  stringstream i_csc_detector_id;
1077  i_csc_detector_id << i_cscdetid;
1078 
1079  stringstream i_ss;
1080  i_ss<<"\t\t hit "<<i_compHit<<" -Muon CSC- detID = "<<i_csc_detector_id.str();
1081 
1082  string i_hitlog = i_ss.str();
1083  i_hitlog = i_hitlog + write_matched_simtracks(i_SimTrackIds);
1084  if (printRtS) edm::LogVerbatim("MuonAssociatorByHitsHelper") << i_hitlog;
1085 
1086  SimTrackIds.insert(SimTrackIds.end(),i_SimTrackIds.begin(),i_SimTrackIds.end());
1087  }
1088  } // if (cscsegment)
1089 
1090  else if (printRtS) edm::LogWarning("MuonAssociatorByHitsHelper")
1091  <<"*** WARNING in MuonAssociatorByHitsHelper::getMatchedIds, CSC TrackingRecHit is neither CSCRecHit2D nor CSCSegment ! ";
1092  }
1093  }
1094 
1095  // RPC Hits
1096  else if (subdet == MuonSubdetId::RPC) {
1097  RPCDetId rpcdetid = RPCDetId(detid);
1098  stringstream rpc_detector_id;
1099  rpc_detector_id << rpcdetid;
1100  if (valid_Hit) hitlog = hitlog+" -Muon RPC- detID = "+rpc_detector_id.str();
1101  else hitlog = hitlog+" *** INVALID ***"+" -Muon RPC- detID = "+rpc_detector_id.str();
1102 
1103  iH++;
1104  SimTrackIds = rpctruth.associateRecHit(*hitp);
1105 
1106  if (valid_Hit) {
1107  n_rpc_valid++;
1108 
1109  if (!SimTrackIds.empty()) {
1110  n_rpc_matched_valid++;
1111  //muon_matchedIds_valid[iH] = SimTrackIds;
1112  muon_matchedIds_valid.push_back (new uint_SimHitIdpr_pair(iH,SimTrackIds));
1113 
1114  }
1115  } else {
1116  n_rpc_INVALID++;
1117 
1118  if (!SimTrackIds.empty()) {
1119  n_rpc_matched_INVALID++;
1120  //muon_matchedIds_INVALID[iH] = SimTrackIds;
1121  muon_matchedIds_INVALID.push_back (new uint_SimHitIdpr_pair(iH,SimTrackIds));
1122  }
1123  }
1124  }
1125 
1126  // GEM Hits
1127  else if (subdet == MuonSubdetId::GEM) {
1128  GEMDetId gemdetid = GEMDetId(detid);
1129  stringstream gem_detector_id;
1130  gem_detector_id << gemdetid;
1131  if (valid_Hit) hitlog = hitlog+" -Muon GEM- detID = "+gem_detector_id.str();
1132  else hitlog = hitlog+" *** INVALID ***"+" -Muon GEM- detID = "+gem_detector_id.str();
1133 
1134  iH++;
1135  SimTrackIds = gemtruth.associateRecHit(*hitp);
1136 
1137  if (valid_Hit) {
1138  n_gem_valid++;
1139 
1140  if (!SimTrackIds.empty()) {
1141  n_gem_matched_valid++;
1142  //muon_matchedIds_valid[iH] = SimTrackIds;
1143  muon_matchedIds_valid.push_back (new uint_SimHitIdpr_pair(iH,SimTrackIds));
1144 
1145  }
1146  } else {
1147  n_gem_INVALID++;
1148 
1149  if (!SimTrackIds.empty()) {
1150  n_gem_matched_INVALID++;
1151  //muon_matchedIds_INVALID[iH] = SimTrackIds;
1152  muon_matchedIds_INVALID.push_back (new uint_SimHitIdpr_pair(iH,SimTrackIds));
1153  }
1154  }
1155 
1156 
1157  } else if (printRtS) edm::LogVerbatim("MuonAssociatorByHitsHelper")
1158  <<"TrackingRecHit "<<iloop<<" *** WARNING *** Unexpected Hit from Detector = "<<det;
1159  }
1160  else continue;
1161 
1162  hitlog = hitlog + write_matched_simtracks(SimTrackIds);
1163 
1164  if (printRtS) edm::LogVerbatim("MuonAssociatorByHitsHelper") << hitlog;
1165  if (printRtS && dumpDT && det==DetId::Muon && subdet==MuonSubdetId::DT) {
1166  edm::LogVerbatim("MuonAssociatorByHitsHelper") <<wireidlog;
1167  for (unsigned int j=0; j<DTSimHits.size(); j++) {
1168  edm::LogVerbatim("MuonAssociatorByHitsHelper") <<DTSimHits[j];
1169  }
1170  }
1171 
1172  } //trackingRecHit loop
1173 }
1174 
1175 int MuonAssociatorByHitsHelper::getShared(MapOfMatchedIds & matchedIds, TrackingParticleCollection::const_iterator trpart) const {
1176  int nshared = 0;
1177 
1178 
1179  // map is indexed over the rechits of the reco::Track (no double-countings allowed)
1180  for (MapOfMatchedIds::const_iterator iRecH=matchedIds.begin(); iRecH!=matchedIds.end(); ++iRecH) {
1181 
1182  // vector of associated simhits associated to the current rechit
1183  std::vector<SimHitIdpr> const & SimTrackIds = (*iRecH).second;
1184 
1185  bool found = false;
1186 
1187  for (std::vector<SimHitIdpr>::const_iterator iSimH=SimTrackIds.begin(); iSimH!=SimTrackIds.end(); ++iSimH) {
1188  uint32_t simtrackId = iSimH->first;
1189  EncodedEventId evtId = iSimH->second;
1190 
1191  // look for shared hits with the given TrackingParticle (looping over component SimTracks)
1192  for (TrackingParticle::g4t_iterator simtrack = trpart->g4Track_begin(); simtrack != trpart->g4Track_end(); ++simtrack) {
1193  if (simtrack->trackId() == simtrackId && simtrack->eventId() == evtId) {
1194  found = true;
1195  break;
1196  }
1197  }
1198 
1199  if (found) {
1200  nshared++;
1201  break;
1202  }
1203  }
1204  }
1205 
1206  return nshared;
1207 }
1208 
1209 std::string MuonAssociatorByHitsHelper::write_matched_simtracks(const std::vector<SimHitIdpr>& SimTrackIds) const {
1210  if (SimTrackIds.empty())
1211  return " *** UNMATCHED ***";
1212 
1213  string hitlog(" matched to SimTrack");
1214 
1215  for(size_t j=0; j<SimTrackIds.size(); j++)
1216  {
1217  char buf[64];
1218  snprintf(buf, 64, " Id:%i/Evt:(%i,%i) ", SimTrackIds[j].first, SimTrackIds[j].second.event(), SimTrackIds[j].second.bunchCrossing());
1219  hitlog += buf;
1220  }
1221  return hitlog;
1222 }
1223 
CSCDetId
CSCDetId
Definition: RecoLocalMuon_RECO.doi:36
TrackingParticleCollection
std::vector< TrackingParticle > TrackingParticleCollection
Definition: HiEvtPlaneFlatProducer.cc:111
edm::ParameterSet::dump
std::string dump(unsigned int indent=0) const
Definition: ParameterSet.cc:924
MuonAssociatorByHitsHelper::write_matched_simtracks
std::string write_matched_simtracks(const std::vector< SimHitIdpr > &) const
Definition: MuonAssociatorByHitsHelper.cc:1209
DTRecSegment4D::phiSegment
const DTChamberRecSegment2D * phiSegment() const
The superPhi segment: 0 if no phi projection available.
Definition: DTRecSegment4D.h:96
MuonAssociatorByHitsHelper::IndexAssociation
std::map< size_t, std::vector< IndexMatch > > IndexAssociation
Definition: MuonAssociatorByHitsHelper.h:54
MuonSubdetId::GEM
static const int GEM
Definition: MuonSubdetId.h:15
MuonAssociatorByHitsHelper::getShared
int getShared(MapOfMatchedIds &matchedIds, TrackingParticleCollection::const_iterator trpart) const
Definition: MuonAssociatorByHitsHelper.cc:1175
mps_fire.resources
tuple resources
Definition: mps_fire.py:58
CSCHitAssociator::associateCSCHitId
std::vector< SimHitIdpr > associateCSCHitId(const CSCRecHit2D *) const
Definition: CSCHitAssociator.cc:35
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:107
TrackerTopology::print
std::string print(DetId detid) const
Definition: TrackerTopology.cc:255
DTHitAssociator::associateHit
std::vector< PSimHit > associateHit(const TrackingRecHit &hit) const
Definition: DTHitAssociator.cc:370
MuonAssociatorByHitsHelper::associateSimToRecoIndices
IndexAssociation associateSimToRecoIndices(const TrackHitsCollection &, const edm::RefVector< TrackingParticleCollection > &, Resources const &) const
Definition: MuonAssociatorByHitsHelper.cc:341
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
DetId::rawId
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:250
RPCHitAssociator::associateRecHit
std::vector< SimHitIdpr > associateRecHit(const TrackingRecHit &hit) const
Definition: RPCHitAssociator.cc:79
MuonAssociatorByHitsHelper::TrackHitsCollection
std::vector< std::pair< trackingRecHit_iterator, trackingRecHit_iterator > > TrackHitsCollection
Definition: MuonAssociatorByHitsHelper.h:34
MuonAssociatorByHitsHelper::MuonAssociatorByHitsHelper
MuonAssociatorByHitsHelper(const edm::ParameterSet &conf)
Definition: MuonAssociatorByHitsHelper.cc:17
MuonSubdetId::CSC
static const int CSC
Definition: MuonSubdetId.h:13
DTRecSegment2D::recHits
virtual std::vector< const TrackingRecHit * > recHits() const
Access to component RecHits (if any)
Definition: DTRecSegment2D.cc:88
j
int j
Definition: DBlmapReader.cc:9
MuonAssociatorByHitsHelper::Resources::diagnostics_
std::function< void(const TrackHitsCollection &, const TrackingParticleCollection &)> diagnostics_
Definition: MuonAssociatorByHitsHelper.h:45
end
#define end
Definition: vmac.h:37
DTRecSegment4D::hasPhi
bool hasPhi() const
Does it have the Phi projection?
Definition: DTRecSegment4D.h:90
DetId::subdetId
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
CSCSegment::recHits
virtual std::vector< const TrackingRecHit * > recHits() const
Access to component RecHits (if any)
Definition: CSCSegment.cc:30
TrackingParticle::g4t_iterator
std::vector< SimTrack >::const_iterator g4t_iterator
Definition: TrackingParticle.h:41
DTRecSegment4D::hasZed
bool hasZed() const
Does it have the Z projection?
Definition: DTRecSegment4D.h:93
DetId
Definition: DetId.h:18
DTRecSegment4D::zSegment
const DTSLRecSegment2D * zSegment() const
The Z segment: 0 if not zed projection available.
Definition: DTRecSegment4D.h:101
DetId::Detector
Detector
Definition: DetId.h:24
TrackingRecHit::isValid
bool isValid() const
Definition: TrackingRecHit.h:139
MuonAssociatorByHitsHelper::uint_SimHitIdpr_pair
std::pair< unsigned int, std::vector< SimHitIdpr > > uint_SimHitIdpr_pair
Definition: MuonAssociatorByHitsHelper.h:32
MuonAssociatorByHitsHelper::associateRecoToSimIndices
IndexAssociation associateRecoToSimIndices(const TrackHitsCollection &, const edm::RefVector< TrackingParticleCollection > &, Resources const &) const
Definition: MuonAssociatorByHitsHelper.cc:59
TrackerHitAssociator::associateHitId
std::vector< SimHitIdpr > associateHitId(const TrackingRecHit &thit) const
Definition: TrackerHitAssociator.cc:309
begin
#define begin
Definition: vmac.h:30
MuonSubdetId::RPC
static const int RPC
Definition: MuonSubdetId.h:14
MuonAssociatorByHitsHelper::MapOfMatchedIds
boost::ptr_vector< uint_SimHitIdpr_pair > MapOfMatchedIds
Definition: MuonAssociatorByHitsHelper.h:33
RPCDetId
RPCDetId
Definition: RecoLocalMuon_RECO.doi:57
edm::RefVector::size
size_type size() const
Size of the RefVector.
Definition: RefVector.h:107
MuonSubdetId::DT
static const int DT
Definition: MuonSubdetId.h:12
TrackingRecHit::geographicalId
DetId geographicalId() const
Definition: TrackingRecHit.h:116
DetId::det
Detector det() const
get the detector field from this detid
Definition: DetId.h:35
GEMHitAssociator::associateRecHit
std::vector< SimHitIdpr > associateRecHit(const TrackingRecHit &hit) const
Definition: GEMHitAssociator.cc:82
DTHitAssociator::associateDTHitId
std::vector< SimHitIdpr > associateDTHitId(const DTRecHit1D *dtrechit) const
Definition: DTHitAssociator.cc:299
trackingRecHit_iterator
TrackingRecHitCollection::base::const_iterator trackingRecHit_iterator
iterator over a vector of reference to TrackingRecHit in the same collection
Definition: TrackingRecHitFwd.h:19
MuonAssociatorByHitsHelper::getMatchedIds
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
Definition: MuonAssociatorByHitsHelper.cc:766
DTRecHit1D::wireId
DTWireId wireId() const
Return the wireId.
Definition: DTRecHit1D.h:107
Generated for CMSSW Reference Manual by doxygen 1.8.5