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