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