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