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