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