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