00001 #include "SimMuon/MCTruth/interface/MuonAssociatorByHits.h"
00002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00003 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00004 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
00005 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00006 #include "DataFormats/DetId/interface/DetId.h"
00007 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00008 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00009 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
00010 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
00011 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
00012 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
00013 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
00014 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
00015 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00016 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00017 #include "DataFormats/DTRecHit/interface/DTRecSegment4D.h"
00018 #include "DataFormats/CSCRecHit/interface/CSCSegment.h"
00019 #include "SimMuon/MCTruth/interface/TrackerMuonHitExtractor.h"
00020 #include <sstream>
00021
00022 using namespace reco;
00023 using namespace std;
00024
00025 MuonAssociatorByHits::MuonAssociatorByHits (const edm::ParameterSet& conf) :
00026 includeZeroHitMuons(conf.getParameter<bool>("includeZeroHitMuons")),
00027 acceptOneStubMatchings(conf.getParameter<bool>("acceptOneStubMatchings")),
00028 UseTracker(conf.getParameter<bool>("UseTracker")),
00029 UseMuon(conf.getParameter<bool>("UseMuon")),
00030 AbsoluteNumberOfHits_track(conf.getParameter<bool>("AbsoluteNumberOfHits_track")),
00031 NHitCut_track(conf.getParameter<unsigned int>("NHitCut_track")),
00032 EfficiencyCut_track(conf.getParameter<double>("EfficiencyCut_track")),
00033 PurityCut_track(conf.getParameter<double>("PurityCut_track")),
00034 AbsoluteNumberOfHits_muon(conf.getParameter<bool>("AbsoluteNumberOfHits_muon")),
00035 NHitCut_muon(conf.getParameter<unsigned int>("NHitCut_muon")),
00036 EfficiencyCut_muon(conf.getParameter<double>("EfficiencyCut_muon")),
00037 PurityCut_muon(conf.getParameter<double>("PurityCut_muon")),
00038 UsePixels(conf.getParameter<bool>("UsePixels")),
00039 UseGrouped(conf.getParameter<bool>("UseGrouped")),
00040 UseSplitting(conf.getParameter<bool>("UseSplitting")),
00041 ThreeHitTracksAreSpecial(conf.getParameter<bool>("ThreeHitTracksAreSpecial")),
00042 dumpDT(conf.getParameter<bool>("dumpDT")),
00043 dumpInputCollections(conf.getParameter<bool>("dumpInputCollections")),
00044 crossingframe(conf.getParameter<bool>("crossingframe")),
00045 simtracksTag(conf.getParameter<edm::InputTag>("simtracksTag")),
00046 simtracksXFTag(conf.getParameter<edm::InputTag>("simtracksXFTag")),
00047 conf_(conf)
00048 {
00049 edm::LogVerbatim("MuonAssociatorByHits") << "constructing MuonAssociatorByHits" << conf_.dump();
00050
00051
00052 if (UseTracker) edm::LogVerbatim("MuonAssociatorByHits")<<"\n UseTracker = TRUE : Tracker SimHits and RecHits WILL be counted";
00053 else edm::LogVerbatim("MuonAssociatorByHits") <<"\n UseTracker = FALSE : Tracker SimHits and RecHits WILL NOT be counted";
00054
00055
00056 if (UseMuon) edm::LogVerbatim("MuonAssociatorByHits")<<" UseMuon = TRUE : Muon SimHits and RecHits WILL be counted";
00057 else edm::LogVerbatim("MuonAssociatorByHits") <<" UseMuon = FALSE : Muon SimHits and RecHits WILL NOT be counted"<<endl;
00058
00059
00060 if (includeZeroHitMuons) {
00061 edm::LogVerbatim("MuonAssociatorByHits")
00062 <<"\n includeZeroHitMuons = TRUE"
00063 <<"\n ==> (re)set NHitCut_muon = 0, PurityCut_muon = 0, EfficiencyCut_muon = 0"<<endl;
00064 NHitCut_muon = 0;
00065 PurityCut_muon = 0.;
00066 EfficiencyCut_muon = 0.;
00067 }
00068
00069 }
00070
00071 MuonAssociatorByHits::~MuonAssociatorByHits()
00072 {
00073 }
00074
00075 RecoToSimCollection
00076 MuonAssociatorByHits::associateRecoToSim( const edm::RefToBaseVector<reco::Track>& tC,
00077 const edm::RefVector<TrackingParticleCollection>& TPCollectionH,
00078 const edm::Event * e, const edm::EventSetup * setup) const{
00079 RecoToSimCollection outputCollection;
00080
00081 TrackHitsCollection tH;
00082 for (edm::RefToBaseVector<reco::Track>::const_iterator it = tC.begin(), ed = tC.end(); it != ed; ++it) {
00083 tH.push_back(std::make_pair((*it)->recHitsBegin(), (*it)->recHitsEnd()));
00084 }
00085
00086 IndexAssociation bareAssoc = associateRecoToSimIndices(tH, TPCollectionH, e, setup);
00087 for (IndexAssociation::const_iterator it = bareAssoc.begin(), ed = bareAssoc.end(); it != ed; ++it) {
00088 for (std::vector<IndexMatch>::const_iterator itma = it->second.begin(), edma = it->second.end(); itma != edma; ++itma) {
00089 outputCollection.insert(tC[it->first], std::make_pair(edm::Ref<TrackingParticleCollection>(TPCollectionH, itma->idx), itma->quality));
00090 }
00091 }
00092
00093 outputCollection.post_insert();
00094 return outputCollection;
00095 }
00096
00097 MuonAssociatorByHits::IndexAssociation
00098 MuonAssociatorByHits::associateRecoToSimIndices(const TrackHitsCollection & tC,
00099 const edm::RefVector<TrackingParticleCollection>& TPCollectionH,
00100 const edm::Event * e, const edm::EventSetup * setup) const{
00101
00102 int tracker_nshared = 0;
00103 int muon_nshared = 0;
00104 int global_nshared = 0;
00105
00106 double tracker_quality = 0;
00107 double tracker_quality_cut;
00108 if (AbsoluteNumberOfHits_track) tracker_quality_cut = static_cast<double>(NHitCut_track);
00109 else tracker_quality_cut = PurityCut_track;
00110
00111 double muon_quality = 0;
00112 double muon_quality_cut;
00113 if (AbsoluteNumberOfHits_muon) muon_quality_cut = static_cast<double>(NHitCut_muon);
00114 else muon_quality_cut = PurityCut_muon;
00115
00116 double global_quality = 0;
00117
00118 MapOfMatchedIds tracker_matchedIds_valid, muon_matchedIds_valid;
00119 MapOfMatchedIds tracker_matchedIds_INVALID, muon_matchedIds_INVALID;
00120
00121 IndexAssociation outputCollection;
00122 bool printRtS(true);
00123
00124
00125 TrackerHitAssociator * trackertruth = new TrackerHitAssociator(*e, conf_);
00126
00127 MuonTruth csctruth(*e,*setup,conf_);
00128
00129 printRtS = true;
00130 DTHitAssociator dttruth(*e,*setup,conf_,printRtS);
00131
00132 RPCHitAssociator rpctruth(*e,*setup,conf_);
00133
00134 TrackingParticleCollection tPC;
00135 if (TPCollectionH.size()!=0) tPC = *(TPCollectionH.product());
00136
00137 if (dumpInputCollections) {
00138
00139 edm::LogVerbatim("MuonAssociatorByHits")<<"\n"<<"reco::Track collection --- size = "<<tC.size();
00140
00141
00142
00143 edm::LogVerbatim("MuonAssociatorByHits")<<"\n"<<"TrackingParticle collection --- size = "<<tPC.size();
00144 int j = 0;
00145 for(TrackingParticleCollection::const_iterator ITER=tPC.begin(); ITER!=tPC.end(); ITER++, j++) {
00146 edm::LogVerbatim("MuonAssociatorByHits")
00147 <<"TrackingParticle "<<j<<", q = "<<ITER->charge()<<", p = "<<ITER->p()
00148 <<", pT = "<<ITER->pt()<<", eta = "<<ITER->eta()<<", phi = "<<ITER->phi();
00149
00150 edm::LogVerbatim("MuonAssociatorByHits")
00151 <<"\t pdg code = "<<ITER->pdgId()<<", made of "<<ITER->trackPSimHit().size()<<" PSimHit"
00152 <<" (in "<<ITER->matchedHit()<<" layers)"
00153 <<" from "<<ITER->g4Tracks().size()<<" SimTrack:";
00154 for (TrackingParticle::g4t_iterator g4T=ITER->g4Track_begin(); g4T!=ITER->g4Track_end(); g4T++) {
00155 edm::LogVerbatim("MuonAssociatorByHits")
00156 <<"\t\t Id:"<<g4T->trackId()<<"/Evt:("<<g4T->eventId().event()<<","<<g4T->eventId().bunchCrossing()<<")";
00157 }
00158 }
00159
00160
00161 edm::Handle<CrossingFrame<SimTrack> > cf_simtracks;
00162 edm::Handle<edm::SimTrackContainer> simTrackCollection;
00163
00164
00165 edm::Handle<CrossingFrame<SimVertex> > cf_simvertices;
00166 edm::Handle<edm::SimVertexContainer> simVertexCollection;
00167
00168 if (crossingframe) {
00169 e->getByLabel(simtracksXFTag,cf_simtracks);
00170 auto_ptr<MixCollection<SimTrack> > SimTk( new MixCollection<SimTrack>(cf_simtracks.product()) );
00171 edm::LogVerbatim("MuonAssociatorByHits")<<"\n"<<"CrossingFrame<SimTrack> collection with InputTag = "<<simtracksXFTag
00172 <<" has size = "<<SimTk->size();
00173 int k = 0;
00174 for (MixCollection<SimTrack>::MixItr ITER=SimTk->begin(); ITER!=SimTk->end(); ITER++, k++) {
00175 edm::LogVerbatim("MuonAssociatorByHits")
00176 <<"SimTrack "<<k
00177 <<" - Id:"<<ITER->trackId()<<"/Evt:("<<ITER->eventId().event()<<","<<ITER->eventId().bunchCrossing()<<")"
00178 <<" pdgId = "<<ITER->type()<<", q = "<<ITER->charge()<<", p = "<<ITER->momentum().P()
00179 <<", pT = "<<ITER->momentum().Pt()<<", eta = "<<ITER->momentum().Eta()<<", phi = "<<ITER->momentum().Phi()
00180 <<"\n * "<<*ITER <<endl;
00181 }
00182 e->getByLabel(simtracksXFTag,cf_simvertices);
00183 auto_ptr<MixCollection<SimVertex> > SimVtx( new MixCollection<SimVertex>(cf_simvertices.product()) );
00184 edm::LogVerbatim("MuonAssociatorByHits")<<"\n"<<"CrossingFrame<SimVertex> collection with InputTag = "<<simtracksXFTag
00185 <<" has size = "<<SimVtx->size();
00186 int kv = 0;
00187 for (MixCollection<SimVertex>::MixItr VITER=SimVtx->begin(); VITER!=SimVtx->end(); VITER++, kv++){
00188 edm::LogVerbatim("MuonAssociatorByHits")
00189 <<"SimVertex "<<kv
00190 << " : "<< *VITER <<endl;
00191 }
00192 }
00193 else {
00194 e->getByLabel(simtracksTag,simTrackCollection);
00195 const edm::SimTrackContainer simTC = *(simTrackCollection.product());
00196 edm::LogVerbatim("MuonAssociatorByHits")<<"\n"<<"SimTrack collection with InputTag = "<<simtracksTag
00197 <<" has size = "<<simTC.size()<<endl;
00198 int k = 0;
00199 for(edm::SimTrackContainer::const_iterator ITER=simTC.begin(); ITER!=simTC.end(); ITER++, k++){
00200 edm::LogVerbatim("MuonAssociatorByHits")
00201 <<"SimTrack "<<k
00202 <<" - Id:"<<ITER->trackId()<<"/Evt:("<<ITER->eventId().event()<<","<<ITER->eventId().bunchCrossing()<<")"
00203 <<" pdgId = "<<ITER->type()<<", q = "<<ITER->charge()<<", p = "<<ITER->momentum().P()
00204 <<", pT = "<<ITER->momentum().Pt()<<", eta = "<<ITER->momentum().Eta()<<", phi = "<<ITER->momentum().Phi()
00205 <<"\n * "<<*ITER <<endl;
00206 }
00207 e->getByLabel("g4SimHits",simVertexCollection);
00208 const edm::SimVertexContainer simVC = *(simVertexCollection.product());
00209 edm::LogVerbatim("MuonAssociatorByHits")<<"\n"<<"SimVertex collection with InputTag = "<<"g4SimHits"
00210 <<" has size = "<<simVC.size()<<endl;
00211 int kv = 0;
00212 for (edm::SimVertexContainer::const_iterator VITER=simVC.begin(); VITER!=simVC.end(); VITER++, kv++){
00213 edm::LogVerbatim("MuonAssociatorByHits")
00214 <<"SimVertex "<<kv
00215 << " : "<< *VITER <<endl;
00216 }
00217 }
00218 }
00219
00220 int tindex=0;
00221 for (TrackHitsCollection::const_iterator track=tC.begin(); track!=tC.end(); track++, tindex++) {
00222 edm::LogVerbatim("MuonAssociatorByHits")
00223 <<"\n"<<"reco::Track "<<tindex
00224 <<", number of RecHits = "<< (track->second - track->first) << "\n";
00225 tracker_matchedIds_valid.clear();
00226 muon_matchedIds_valid.clear();
00227
00228 tracker_matchedIds_INVALID.clear();
00229 muon_matchedIds_INVALID.clear();
00230
00231 bool this_track_matched = false;
00232 int n_matching_simhits = 0;
00233
00234
00235 int n_all = 0;
00236 int n_tracker_all = 0;
00237 int n_muon_all = 0;
00238 int n_dt_all = 0;
00239 int n_csc_all = 0;
00240 int n_rpc_all = 0;
00241
00242 int n_valid = 0;
00243 int n_tracker_valid = 0;
00244 int n_muon_valid = 0;
00245 int n_dt_valid = 0;
00246 int n_csc_valid = 0;
00247 int n_rpc_valid = 0;
00248
00249 int n_matched_valid = 0;
00250 int n_tracker_matched_valid = 0;
00251 int n_muon_matched_valid = 0;
00252 int n_dt_matched_valid = 0;
00253 int n_csc_matched_valid = 0;
00254 int n_rpc_matched_valid = 0;
00255
00256 int n_INVALID = 0;
00257 int n_tracker_INVALID = 0;
00258 int n_muon_INVALID = 0;
00259 int n_dt_INVALID = 0;
00260 int n_csc_INVALID = 0;
00261 int n_rpc_INVALID = 0;
00262
00263 int n_matched_INVALID = 0;
00264 int n_tracker_matched_INVALID = 0;
00265 int n_muon_matched_INVALID = 0;
00266 int n_dt_matched_INVALID = 0;
00267 int n_csc_matched_INVALID = 0;
00268 int n_rpc_matched_INVALID = 0;
00269
00270 printRtS = true;
00271 getMatchedIds(tracker_matchedIds_valid, muon_matchedIds_valid,
00272 tracker_matchedIds_INVALID, muon_matchedIds_INVALID,
00273 n_tracker_valid, n_dt_valid, n_csc_valid, n_rpc_valid,
00274 n_tracker_matched_valid, n_dt_matched_valid, n_csc_matched_valid, n_rpc_matched_valid,
00275 n_tracker_INVALID, n_dt_INVALID, n_csc_INVALID, n_rpc_INVALID,
00276 n_tracker_matched_INVALID, n_dt_matched_INVALID, n_csc_matched_INVALID, n_rpc_matched_INVALID,
00277 track->first, track->second,
00278 trackertruth, dttruth, csctruth, rpctruth,
00279 printRtS);
00280
00281 n_matching_simhits = tracker_matchedIds_valid.size() + muon_matchedIds_valid.size() +
00282 tracker_matchedIds_INVALID.size() +muon_matchedIds_INVALID.size();
00283
00284 n_muon_valid = n_dt_valid + n_csc_valid + n_rpc_valid;
00285 n_valid = n_tracker_valid + n_muon_valid;
00286 n_muon_INVALID = n_dt_INVALID + n_csc_INVALID + n_rpc_INVALID;
00287 n_INVALID = n_tracker_INVALID + n_muon_INVALID;
00288
00289
00290 n_tracker_all = n_tracker_valid + n_tracker_INVALID;
00291 n_dt_all = n_dt_valid + n_dt_INVALID;
00292 n_csc_all = n_csc_valid + n_csc_INVALID;
00293 n_rpc_all = n_rpc_valid + n_rpc_INVALID;
00294 n_muon_all = n_muon_valid + n_muon_INVALID;
00295 n_all = n_valid + n_INVALID;
00296
00297 n_muon_matched_valid = n_dt_matched_valid + n_csc_matched_valid + n_rpc_matched_valid;
00298 n_matched_valid = n_tracker_matched_valid + n_muon_matched_valid;
00299 n_muon_matched_INVALID = n_dt_matched_INVALID + n_csc_matched_INVALID + n_rpc_matched_INVALID;
00300 n_matched_INVALID = n_tracker_matched_INVALID + n_muon_matched_INVALID;
00301
00302
00303 int n_tracker_selected_hits = n_tracker_valid;
00304 int n_muon_selected_hits = n_muon_valid;
00305 int n_dt_selected_hits = n_dt_valid;
00306 int n_csc_selected_hits = n_csc_valid;
00307 int n_rpc_selected_hits = n_rpc_valid;
00308
00309
00310 int n_tracker_matched = n_tracker_matched_valid;
00311 int n_muon_matched = n_muon_matched_valid;
00312 int n_dt_matched = n_dt_matched_valid;
00313 int n_csc_matched = n_csc_matched_valid;
00314 int n_rpc_matched = n_rpc_matched_valid;
00315
00316 std::string InvMuonHits, ZeroHitMuon;
00317
00318 if (includeZeroHitMuons && n_muon_valid==0 && n_muon_INVALID!=0) {
00319
00320
00321 InvMuonHits = " ***INVALID MUON HITS***";
00322 ZeroHitMuon = " ***ZERO-HIT MUON***";
00323
00324 n_muon_selected_hits = n_muon_INVALID;
00325 n_dt_selected_hits = n_dt_INVALID;
00326 n_csc_selected_hits = n_csc_INVALID;
00327 n_rpc_selected_hits = n_rpc_INVALID;
00328
00329 n_muon_matched = n_muon_matched_INVALID;
00330 n_dt_matched = n_dt_matched_INVALID;
00331 n_csc_matched = n_csc_matched_INVALID;
00332 n_rpc_matched = n_rpc_matched_INVALID;
00333 }
00334
00335 int n_selected_hits = n_tracker_selected_hits + n_muon_selected_hits;
00336 int n_matched = n_tracker_matched + n_muon_matched;
00337
00338 edm::LogVerbatim("MuonAssociatorByHits")
00339 <<"\n"<<"# TrackingRecHits: "<<(track->second - track->first)
00340 <<"\n"<< "# used RecHits = " << n_all <<" ("<<n_tracker_all<<"/"
00341 <<n_dt_all<<"/"<<n_csc_all<<"/"<<n_rpc_all<<" in Tracker/DT/CSC/RPC)"<<", obtained from " << n_matching_simhits << " SimHits"
00342 <<"\n"<< "# selected RecHits = " <<n_selected_hits <<" (" <<n_tracker_selected_hits<<"/"
00343 <<n_dt_selected_hits<<"/"<<n_csc_selected_hits<<"/"<<n_rpc_selected_hits<<" in Tracker/DT/CSC/RPC)"<<InvMuonHits
00344 <<"\n"<< "# matched RecHits = " <<n_matched<<" ("<<n_tracker_matched<<"/"
00345 <<n_dt_matched<<"/"<<n_csc_matched<<"/"<<n_rpc_matched<<" in Tracker/DT/CSC/RPC)";
00346
00347 if (n_all>0 && n_matching_simhits == 0)
00348 edm::LogWarning("MuonAssociatorByHits")
00349 <<"*** WARNING in MuonAssociatorByHits::associateRecoToSim: no matching PSimHit found for this reco::Track !";
00350
00351 if (n_matching_simhits != 0) {
00352 edm::LogVerbatim("MuonAssociatorByHits")
00353 <<"\n"<< "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
00354 <<"\n"<< "reco::Track "<<tindex<<ZeroHitMuon
00355 <<"\n\t"<< "made of "<<n_selected_hits<<" selected RecHits (tracker:"<<n_tracker_selected_hits<<"/muons:"<<n_muon_selected_hits<<")";
00356
00357 int tpindex = 0;
00358 for (TrackingParticleCollection::const_iterator trpart = tPC.begin(); trpart != tPC.end(); ++trpart, ++tpindex) {
00359 tracker_nshared = getShared(tracker_matchedIds_valid, trpart);
00360 muon_nshared = getShared(muon_matchedIds_valid, trpart);
00361
00362 if (includeZeroHitMuons && n_muon_valid==0 && n_muon_INVALID!=0)
00363 muon_nshared = getShared(muon_matchedIds_INVALID, trpart);
00364
00365 global_nshared = tracker_nshared + muon_nshared;
00366
00367 if (AbsoluteNumberOfHits_track) tracker_quality = static_cast<double>(tracker_nshared);
00368 else if(n_tracker_selected_hits != 0) tracker_quality = (static_cast<double>(tracker_nshared)/static_cast<double>(n_tracker_selected_hits));
00369 else tracker_quality = 0;
00370
00371 if (AbsoluteNumberOfHits_muon) muon_quality = static_cast<double>(muon_nshared);
00372 else if(n_muon_selected_hits != 0) muon_quality = (static_cast<double>(muon_nshared)/static_cast<double>(n_muon_selected_hits));
00373 else muon_quality = 0;
00374
00375
00376 if (n_selected_hits != 0) {
00377 if (AbsoluteNumberOfHits_muon && AbsoluteNumberOfHits_track)
00378 global_quality = global_nshared;
00379 else
00380 global_quality = (static_cast<double>(global_nshared)/static_cast<double>(n_selected_hits));
00381 } else global_quality = 0;
00382
00383 bool trackerOk = false;
00384 if (n_tracker_selected_hits != 0) {
00385 if (tracker_quality > tracker_quality_cut) trackerOk = true;
00386
00387 if (ThreeHitTracksAreSpecial && n_tracker_selected_hits==3 && tracker_nshared<3) trackerOk = false;
00388 }
00389
00390 bool muonOk = false;
00391 if (n_muon_selected_hits != 0) {
00392 if (muon_quality > muon_quality_cut) muonOk = true;
00393 }
00394
00395
00396 bool matchOk = trackerOk || muonOk;
00397
00398
00399 if (!acceptOneStubMatchings && n_tracker_selected_hits!=0 && n_muon_selected_hits!=0)
00400 matchOk = trackerOk && muonOk;
00401
00402 if (matchOk) {
00403
00404 outputCollection[tindex].push_back(IndexMatch(tpindex, global_quality));
00405 this_track_matched = true;
00406
00407 edm::LogVerbatim("MuonAssociatorByHits")
00408 << "\n\t"<<" **MATCHED** with quality = "<<global_quality<<" (tracker: "<<tracker_quality<<" / muon: "<<muon_quality<<")"
00409 << "\n\t"<<" N shared hits = "<<global_nshared<<" (tracker: "<<tracker_nshared<<" / muon: "<<muon_nshared<<")"
00410 <<"\n"<< " to: TrackingParticle " <<tpindex<<", q = "<<(*trpart).charge()<<", p = "<<(*trpart).p()
00411 <<", pT = "<<(*trpart).pt()<<", eta = "<<(*trpart).eta()<<", phi = "<<(*trpart).phi()
00412 <<"\n\t"<< " pdg code = "<<(*trpart).pdgId()<<", made of "<<(*trpart).trackPSimHit().size()<<" PSimHits"
00413 <<" from "<<(*trpart).g4Tracks().size()<<" SimTrack:";
00414 for(TrackingParticle::g4t_iterator g4T=(*trpart).g4Track_begin();
00415 g4T!=(*trpart).g4Track_end();
00416 ++g4T) {
00417 edm::LogVerbatim("MuonAssociatorByHits")
00418 <<"\t"<< " Id:"<<(*g4T).trackId()<<"/Evt:("<<(*g4T).eventId().event()<<","<<(*g4T).eventId().bunchCrossing()<<")";
00419 }
00420 }
00421 else {
00422
00423 if (global_nshared != 0)
00424 edm::LogVerbatim("MuonAssociatorByHits")
00425 <<"\n\t"<<" NOT matched to TrackingParticle "<<tpindex
00426 << " with quality = "<<global_quality<<" (tracker: "<<tracker_quality<<" / muon: "<<muon_quality<<")"
00427 << "\n"<< " N shared hits = "<<global_nshared<<" (tracker: "<<tracker_nshared<<" / muon: "<<muon_nshared<<")";
00428 }
00429
00430 }
00431
00432 if (!this_track_matched) {
00433 edm::LogVerbatim("MuonAssociatorByHits")
00434 <<"\n"<<" NOT matched to any TrackingParticle";
00435 }
00436
00437 edm::LogVerbatim("MuonAssociatorByHits")
00438 <<"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"<<"\n";
00439
00440 }
00441
00442 }
00443
00444 if (!tC.size())
00445 edm::LogVerbatim("MuonAssociatorByHits")<<"0 reconstructed tracks (-->> 0 associated !)";
00446
00447 delete trackertruth;
00448 for (IndexAssociation::iterator it = outputCollection.begin(), ed = outputCollection.end(); it != ed; ++it) {
00449 std::sort(it->second.begin(), it->second.end());
00450 }
00451 return outputCollection;
00452 }
00453
00454
00455 SimToRecoCollection
00456 MuonAssociatorByHits::associateSimToReco( const edm::RefToBaseVector<reco::Track>& tC,
00457 const edm::RefVector<TrackingParticleCollection>& TPCollectionH,
00458 const edm::Event * e, const edm::EventSetup * setup) const{
00459
00460 SimToRecoCollection outputCollection;
00461 TrackHitsCollection tH;
00462 for (edm::RefToBaseVector<reco::Track>::const_iterator it = tC.begin(), ed = tC.end(); it != ed; ++it) {
00463 tH.push_back(std::make_pair((*it)->recHitsBegin(), (*it)->recHitsEnd()));
00464 }
00465
00466 IndexAssociation bareAssoc = associateSimToRecoIndices(tH, TPCollectionH, e, setup);
00467 for (IndexAssociation::const_iterator it = bareAssoc.begin(), ed = bareAssoc.end(); it != ed; ++it) {
00468 for (std::vector<IndexMatch>::const_iterator itma = it->second.begin(), edma = it->second.end(); itma != edma; ++itma) {
00469 outputCollection.insert(edm::Ref<TrackingParticleCollection>(TPCollectionH, it->first),
00470 std::make_pair(tC[itma->idx], itma->quality));
00471 }
00472 }
00473
00474 outputCollection.post_insert();
00475 return outputCollection;
00476 }
00477
00478 MuonAssociatorByHits::IndexAssociation
00479 MuonAssociatorByHits::associateSimToRecoIndices( const TrackHitsCollection & tC,
00480 const edm::RefVector<TrackingParticleCollection>& TPCollectionH,
00481 const edm::Event * e, const edm::EventSetup * setup) const{
00482
00483
00484 int tracker_nshared = 0;
00485 int muon_nshared = 0;
00486 int global_nshared = 0;
00487
00488 double tracker_quality = 0;
00489 double tracker_quality_cut;
00490 if (AbsoluteNumberOfHits_track) tracker_quality_cut = static_cast<double>(NHitCut_track);
00491 else tracker_quality_cut = EfficiencyCut_track;
00492
00493 double muon_quality = 0;
00494 double muon_quality_cut;
00495 if (AbsoluteNumberOfHits_muon) muon_quality_cut = static_cast<double>(NHitCut_muon);
00496 else muon_quality_cut = EfficiencyCut_muon;
00497
00498 double global_quality = 0;
00499
00500 double tracker_purity = 0;
00501 double muon_purity = 0;
00502 double global_purity = 0;
00503
00504 MapOfMatchedIds tracker_matchedIds_valid, muon_matchedIds_valid;
00505 MapOfMatchedIds tracker_matchedIds_INVALID, muon_matchedIds_INVALID;
00506
00507 IndexAssociation outputCollection;
00508
00509 bool printRtS(true);
00510
00511
00512 TrackerHitAssociator * trackertruth = new TrackerHitAssociator(*e, conf_);
00513
00514 MuonTruth csctruth(*e,*setup,conf_);
00515
00516 printRtS = false;
00517 DTHitAssociator dttruth(*e,*setup,conf_,printRtS);
00518
00519 RPCHitAssociator rpctruth(*e,*setup,conf_);
00520
00521 TrackingParticleCollection tPC;
00522 if (TPCollectionH.size()!=0) tPC = *(TPCollectionH.product());
00523
00524 bool any_trackingParticle_matched = false;
00525
00526 int tindex=0;
00527 for (TrackHitsCollection::const_iterator track=tC.begin(); track!=tC.end(); track++, tindex++) {
00528 if (printRtS) edm::LogVerbatim("MuonAssociatorByHits")
00529 <<"\n"<<"reco::Track "<<tindex
00530 <<", number of RecHits = "<< (track->second - track->first) << "\n";
00531
00532 tracker_matchedIds_valid.clear();
00533 muon_matchedIds_valid.clear();
00534
00535 tracker_matchedIds_INVALID.clear();
00536 muon_matchedIds_INVALID.clear();
00537
00538 int n_matching_simhits = 0;
00539
00540
00541 int n_all = 0;
00542 int n_tracker_all = 0;
00543 int n_muon_all = 0;
00544 int n_dt_all = 0;
00545 int n_csc_all = 0;
00546 int n_rpc_all = 0;
00547
00548 int n_valid = 0;
00549 int n_tracker_valid = 0;
00550 int n_muon_valid = 0;
00551 int n_dt_valid = 0;
00552 int n_csc_valid = 0;
00553 int n_rpc_valid = 0;
00554
00555 int n_matched_valid = 0;
00556 int n_tracker_matched_valid = 0;
00557 int n_muon_matched_valid = 0;
00558 int n_dt_matched_valid = 0;
00559 int n_csc_matched_valid = 0;
00560 int n_rpc_matched_valid = 0;
00561
00562 int n_INVALID = 0;
00563 int n_tracker_INVALID = 0;
00564 int n_muon_INVALID = 0;
00565 int n_dt_INVALID = 0;
00566 int n_csc_INVALID = 0;
00567 int n_rpc_INVALID = 0;
00568
00569 int n_matched_INVALID = 0;
00570 int n_tracker_matched_INVALID = 0;
00571 int n_muon_matched_INVALID = 0;
00572 int n_dt_matched_INVALID = 0;
00573 int n_csc_matched_INVALID = 0;
00574 int n_rpc_matched_INVALID = 0;
00575
00576 printRtS = false;
00577 getMatchedIds(tracker_matchedIds_valid, muon_matchedIds_valid,
00578 tracker_matchedIds_INVALID, muon_matchedIds_INVALID,
00579 n_tracker_valid, n_dt_valid, n_csc_valid, n_rpc_valid,
00580 n_tracker_matched_valid, n_dt_matched_valid, n_csc_matched_valid, n_rpc_matched_valid,
00581 n_tracker_INVALID, n_dt_INVALID, n_csc_INVALID, n_rpc_INVALID,
00582 n_tracker_matched_INVALID, n_dt_matched_INVALID, n_csc_matched_INVALID, n_rpc_matched_INVALID,
00583 track->first, track->second,
00584 trackertruth, dttruth, csctruth, rpctruth,
00585 printRtS);
00586
00587 n_matching_simhits = tracker_matchedIds_valid.size() + muon_matchedIds_valid.size() +
00588 tracker_matchedIds_INVALID.size() +muon_matchedIds_INVALID.size();
00589
00590 n_muon_valid = n_dt_valid + n_csc_valid + n_rpc_valid;
00591 n_valid = n_tracker_valid + n_muon_valid;
00592 n_muon_INVALID = n_dt_INVALID + n_csc_INVALID + n_rpc_INVALID;
00593 n_INVALID = n_tracker_INVALID + n_muon_INVALID;
00594
00595
00596 n_tracker_all = n_tracker_valid + n_tracker_INVALID;
00597 n_dt_all = n_dt_valid + n_dt_INVALID;
00598 n_csc_all = n_csc_valid + n_csc_INVALID;
00599 n_rpc_all = n_rpc_valid + n_rpc_INVALID;
00600 n_muon_all = n_muon_valid + n_muon_INVALID;
00601 n_all = n_valid + n_INVALID;
00602
00603 n_muon_matched_valid = n_dt_matched_valid + n_csc_matched_valid + n_rpc_matched_valid;
00604 n_matched_valid = n_tracker_matched_valid + n_muon_matched_valid;
00605 n_muon_matched_INVALID = n_dt_matched_INVALID + n_csc_matched_INVALID + n_rpc_matched_INVALID;
00606 n_matched_INVALID = n_tracker_matched_INVALID + n_muon_matched_INVALID;
00607
00608
00609 int n_tracker_selected_hits = n_tracker_valid;
00610 int n_muon_selected_hits = n_muon_valid;
00611 int n_dt_selected_hits = n_dt_valid;
00612 int n_csc_selected_hits = n_csc_valid;
00613 int n_rpc_selected_hits = n_rpc_valid;
00614
00615
00616 int n_tracker_matched = n_tracker_matched_valid;
00617 int n_muon_matched = n_muon_matched_valid;
00618 int n_dt_matched = n_dt_matched_valid;
00619 int n_csc_matched = n_csc_matched_valid;
00620 int n_rpc_matched = n_rpc_matched_valid;
00621
00622 std::string InvMuonHits, ZeroHitMuon;
00623
00624 if (includeZeroHitMuons && n_muon_valid==0 && n_muon_INVALID!=0) {
00625
00626
00627 InvMuonHits = " ***INVALID MUON HITS***";
00628 ZeroHitMuon = " ***ZERO-HIT MUON***";
00629
00630 n_muon_selected_hits = n_muon_INVALID;
00631 n_dt_selected_hits = n_dt_INVALID;
00632 n_csc_selected_hits = n_csc_INVALID;
00633 n_rpc_selected_hits = n_rpc_INVALID;
00634
00635 n_muon_matched = n_muon_matched_INVALID;
00636 n_dt_matched = n_dt_matched_INVALID;
00637 n_csc_matched = n_csc_matched_INVALID;
00638 n_rpc_matched = n_rpc_matched_INVALID;
00639 }
00640
00641 int n_selected_hits = n_tracker_selected_hits + n_muon_selected_hits;
00642 int n_matched = n_tracker_matched + n_muon_matched;
00643
00644 if (printRtS) edm::LogVerbatim("MuonAssociatorByHits")
00645 <<"\n"<<"# TrackingRecHits: "<<(track->second - track->first)
00646 <<"\n"<< "# used RecHits = " <<n_all <<" ("<<n_tracker_all<<"/"
00647 <<n_dt_all<<"/"<<n_csc_all<<"/"<<n_rpc_all<<" in Tracker/DT/CSC/RPC)"<<", obtained from " << n_matching_simhits << " SimHits"
00648 <<"\n"<< "# selected RecHits = " <<n_selected_hits <<" (" <<n_tracker_selected_hits<<"/"
00649 <<n_dt_selected_hits<<"/"<<n_csc_selected_hits<<"/"<<n_rpc_selected_hits<<" in Tracker/DT/CSC/RPC)"<<InvMuonHits
00650 <<"\n"<< "# matched RecHits = " <<n_matched<<" ("<<n_tracker_matched<<"/"
00651 <<n_dt_matched<<"/"<<n_csc_matched<<"/"<<n_rpc_matched<<" in Tracker/DT/CSC/RPC)";
00652
00653 if (printRtS && n_all>0 && n_matching_simhits==0)
00654 edm::LogWarning("MuonAssociatorByHits")
00655 <<"*** WARNING in MuonAssociatorByHits::associateSimToReco: no matching PSimHit found for this reco::Track !";
00656
00657 if (n_matching_simhits != 0) {
00658 int tpindex =0;
00659 for (TrackingParticleCollection::const_iterator trpart = tPC.begin(); trpart != tPC.end(); ++trpart, ++tpindex) {
00660
00661 int n_tracker_simhits = 0;
00662 int n_tracker_recounted_simhits = 0;
00663 int n_muon_simhits = 0;
00664 int n_global_simhits = 0;
00665 std::vector<PSimHit> tphits;
00666
00667 int n_tracker_selected_simhits = 0;
00668 int n_muon_selected_simhits = 0;
00669 int n_global_selected_simhits = 0;
00670
00671
00672 tracker_nshared = getShared(tracker_matchedIds_valid, trpart);
00673 muon_nshared = getShared(muon_matchedIds_valid, trpart);
00674
00675 if (includeZeroHitMuons && n_muon_valid==0 && n_muon_INVALID!=0)
00676 muon_nshared = getShared(muon_matchedIds_INVALID, trpart);
00677
00678 global_nshared = tracker_nshared + muon_nshared;
00679 if (global_nshared == 0) continue;
00680
00681 for(std::vector<PSimHit>::const_iterator TPhit = trpart->pSimHit_begin(); TPhit != trpart->pSimHit_end(); TPhit++) {
00682 DetId dId = DetId(TPhit->detUnitId());
00683 DetId::Detector detector = dId.det();
00684
00685 if (detector == DetId::Tracker) {
00686 n_tracker_simhits++;
00687
00688 unsigned int subdetId = static_cast<unsigned int>(dId.subdetId());
00689 if (!UsePixels && (subdetId==PixelSubdetector::PixelBarrel || subdetId==PixelSubdetector::PixelEndcap) )
00690 continue;
00691
00692 SiStripDetId* stripDetId = 0;
00693 if (subdetId==SiStripDetId::TIB||subdetId==SiStripDetId::TOB||
00694 subdetId==SiStripDetId::TID||subdetId==SiStripDetId::TEC)
00695 stripDetId= new SiStripDetId(dId);
00696
00697 bool newhit = true;
00698 for(std::vector<PSimHit>::const_iterator TPhitOK = tphits.begin(); TPhitOK != tphits.end(); TPhitOK++) {
00699 DetId dIdOK = DetId(TPhitOK->detUnitId());
00700
00701 if (!UseGrouped && !UseSplitting)
00702 if (LayerFromDetid(dId)==LayerFromDetid(dIdOK) &&
00703 dId.subdetId()==dIdOK.subdetId()) newhit = false;
00704
00705 if (!UseGrouped && UseSplitting)
00706 if (LayerFromDetid(dId)==LayerFromDetid(dIdOK) &&
00707 dId.subdetId()==dIdOK.subdetId() &&
00708 (stripDetId==0 || stripDetId->partnerDetId()!=dIdOK.rawId()))
00709 newhit = false;
00710
00711 if (UseGrouped && !UseSplitting)
00712 if (LayerFromDetid(dId)==LayerFromDetid(dIdOK) &&
00713 dId.subdetId()==dIdOK.subdetId() &&
00714 stripDetId!=0 && stripDetId->partnerDetId()==dIdOK.rawId())
00715 newhit = false;
00716
00717 if (UseGrouped && UseSplitting)
00718 newhit = true;
00719 }
00720 if (newhit) {
00721 tphits.push_back(*TPhit);
00722 }
00723 delete stripDetId;
00724 }
00725 else if (detector == DetId::Muon) {
00726 n_muon_simhits++;
00727
00728
00729 if (dId.subdetId() == MuonSubdetId::CSC) {
00730 if (csctruth.cscBadChambers->isInBadChamber(CSCDetId(dId))) {
00731
00732 n_muon_simhits--;
00733 }
00734 }
00735
00736 }
00737 }
00738
00739 n_tracker_recounted_simhits = tphits.size();
00740
00741 if (trpart->trackPSimHit().empty()) {
00742
00743 n_tracker_recounted_simhits = tracker_nshared;
00744 n_muon_simhits = muon_nshared;
00745 }
00746 n_global_simhits = n_tracker_recounted_simhits + n_muon_simhits;
00747
00748 if (UseMuon) {
00749 n_muon_selected_simhits = n_muon_simhits;
00750 n_global_selected_simhits = n_muon_selected_simhits;
00751 }
00752 if (UseTracker) {
00753 n_tracker_selected_simhits = n_tracker_recounted_simhits;
00754 n_global_selected_simhits += n_tracker_selected_simhits;
00755 }
00756
00757 if (AbsoluteNumberOfHits_track) tracker_quality = static_cast<double>(tracker_nshared);
00758 else if (n_tracker_selected_simhits!=0)
00759 tracker_quality = static_cast<double>(tracker_nshared)/static_cast<double>(n_tracker_selected_simhits);
00760 else tracker_quality = 0;
00761
00762 if (AbsoluteNumberOfHits_muon) muon_quality = static_cast<double>(muon_nshared);
00763 else if (n_muon_selected_simhits!=0)
00764 muon_quality = static_cast<double>(muon_nshared)/static_cast<double>(n_muon_selected_simhits);
00765 else muon_quality = 0;
00766
00767
00768 if (n_global_selected_simhits != 0) {
00769 if (AbsoluteNumberOfHits_muon && AbsoluteNumberOfHits_track)
00770 global_quality = global_nshared;
00771 else
00772 global_quality = static_cast<double>(global_nshared)/static_cast<double>(n_global_selected_simhits);
00773 }
00774 else global_quality = 0;
00775
00776
00777 if (n_selected_hits != 0) {
00778 if (AbsoluteNumberOfHits_muon && AbsoluteNumberOfHits_track)
00779 global_purity = global_nshared;
00780 else
00781 global_purity = static_cast<double>(global_nshared)/static_cast<double>(n_selected_hits);
00782 }
00783 else global_purity = 0;
00784
00785 bool trackerOk = false;
00786 if (n_tracker_selected_hits != 0) {
00787 if (tracker_quality > tracker_quality_cut) trackerOk = true;
00788
00789 tracker_purity = static_cast<double>(tracker_nshared)/static_cast<double>(n_tracker_selected_hits);
00790 if (AbsoluteNumberOfHits_track) tracker_purity = static_cast<double>(tracker_nshared);
00791
00792 if ((!AbsoluteNumberOfHits_track) && tracker_purity <= PurityCut_track) trackerOk = false;
00793
00794
00795 if (ThreeHitTracksAreSpecial && n_tracker_selected_hits==3 && tracker_nshared<3) trackerOk = false;
00796 }
00797
00798 bool muonOk = false;
00799 if (n_muon_selected_hits != 0) {
00800 if (muon_quality > muon_quality_cut) muonOk = true;
00801
00802 muon_purity = static_cast<double>(muon_nshared)/static_cast<double>(n_muon_selected_hits);
00803 if (AbsoluteNumberOfHits_muon) muon_purity = static_cast<double>(muon_nshared);
00804
00805 if ((!AbsoluteNumberOfHits_muon) && muon_purity <= PurityCut_muon) muonOk = false;
00806 }
00807
00808
00809 bool matchOk = trackerOk || muonOk;
00810
00811
00812 if (!acceptOneStubMatchings && n_tracker_selected_hits!=0 && n_muon_selected_hits!=0)
00813 matchOk = trackerOk && muonOk;
00814
00815 if (matchOk) {
00816
00817 outputCollection[tpindex].push_back(IndexMatch(tindex,global_quality));
00818 any_trackingParticle_matched = true;
00819
00820 edm::LogVerbatim("MuonAssociatorByHits")
00821 <<"************************************************************************************************************************"
00822 <<"\n"<< "TrackingParticle " << tpindex <<", q = "<<(*trpart).charge()<<", p = "<<(*trpart).p()
00823 <<", pT = "<<(*trpart).pt()<<", eta = "<<(*trpart).eta()<<", phi = "<<(*trpart).phi()
00824 <<"\n"<<" pdg code = "<<(*trpart).pdgId()
00825 <<", made of "<<(*trpart).trackPSimHit().size()<<" PSimHits, recounted "<<n_global_simhits<<" PSimHits"
00826 <<" (tracker:"<<n_tracker_recounted_simhits<<"/muons:"<<n_muon_simhits<<")"
00827 <<", from "<<(*trpart).g4Tracks().size()<<" SimTrack:";
00828 for(TrackingParticle::g4t_iterator g4T=(*trpart).g4Track_begin();
00829 g4T!=(*trpart).g4Track_end();
00830 ++g4T) {
00831 edm::LogVerbatim("MuonAssociatorByHits")
00832 <<" Id:"<<(*g4T).trackId()<<"/Evt:("<<(*g4T).eventId().event()<<","<<(*g4T).eventId().bunchCrossing()<<")";
00833 }
00834 edm::LogVerbatim("MuonAssociatorByHits")
00835 <<"\t selected "<<n_global_selected_simhits<<" PSimHits"
00836 <<" (tracker:"<<n_tracker_selected_simhits<<"/muons:"<<n_muon_selected_simhits<<")"
00837 << "\n\t **MATCHED** with quality = "<<global_quality<<" (tracker: "<<tracker_quality<<" / muon: "<<muon_quality<<")"
00838 << "\n\t and purity = "<<global_purity<<" (tracker: "<<tracker_purity<<" / muon: "<<muon_purity<<")"
00839 << "\n\t N shared hits = "<<global_nshared<<" (tracker: "<<tracker_nshared<<" / muon: "<<muon_nshared<<")"
00840 <<"\n" <<" to: reco::Track "<<tindex<<ZeroHitMuon
00841 <<"\n\t"<< " made of "<<n_selected_hits<<" RecHits (tracker:"<<n_tracker_valid<<"/muons:"<<n_muon_selected_hits<<")";
00842 }
00843 else {
00844
00845 if (global_nshared != 0) {
00846 if (printRtS) edm::LogVerbatim("MuonAssociatorByHits")
00847 <<"************************************************************************************************************************"
00848 <<"\n"<<"TrackingParticle " << tpindex <<", q = "<<(*trpart).charge()<<", p = "<<(*trpart).p()
00849 <<", pT = "<<(*trpart).pt()<<", eta = "<<(*trpart).eta()<<", phi = "<<(*trpart).phi()
00850 <<"\n"<<" pdg code = "<<(*trpart).pdgId()
00851 <<", made of "<<(*trpart).trackPSimHit().size()<<" PSimHits, recounted "<<n_global_simhits<<" PSimHits"
00852 <<" (tracker:"<<n_tracker_recounted_simhits<<"/muons:"<<n_muon_simhits<<")"
00853 <<", from "<<(*trpart).g4Tracks().size()<<" SimTrack:";
00854 for(TrackingParticle::g4t_iterator g4T=(*trpart).g4Track_begin();
00855 g4T!=(*trpart).g4Track_end();
00856 ++g4T) {
00857 if (printRtS) edm::LogVerbatim("MuonAssociatorByHits")
00858 <<" Id:"<<(*g4T).trackId()<<"/Evt:("<<(*g4T).eventId().event()<<","<<(*g4T).eventId().bunchCrossing()<<")";
00859 }
00860 if (printRtS) edm::LogVerbatim("MuonAssociatorByHits")
00861 <<"\t selected "<<n_global_selected_simhits<<" PSimHits"
00862 <<" (tracker:"<<n_tracker_selected_simhits<<"/muons:"<<n_muon_selected_simhits<<")"
00863 <<"\n\t NOT matched to reco::Track "<<tindex<<ZeroHitMuon
00864 <<" with quality = "<<global_quality<<" (tracker: "<<tracker_quality<<" / muon: "<<muon_quality<<")"
00865 <<"\n\t and purity = "<<global_purity<<" (tracker: "<<tracker_purity<<" / muon: "<<muon_purity<<")"
00866 <<"\n\t N shared hits = "<<global_nshared<<" (tracker: "<<tracker_nshared<<" / muon: "<<muon_nshared<<")";
00867 }
00868 }
00869 }
00870 }
00871 }
00872
00873 if (!any_trackingParticle_matched) {
00874 edm::LogVerbatim("MuonAssociatorByHits")
00875 <<"\n"
00876 <<"************************************************************************************************************************"
00877 << "\n NO TrackingParticle associated to ANY input reco::Track ! \n"
00878 <<"************************************************************************************************************************"<<"\n";
00879 } else {
00880 edm::LogVerbatim("MuonAssociatorByHits")
00881 <<"************************************************************************************************************************"<<"\n";
00882 }
00883
00884 delete trackertruth;
00885 for (IndexAssociation::iterator it = outputCollection.begin(), ed = outputCollection.end(); it != ed; ++it) {
00886 std::sort(it->second.begin(), it->second.end());
00887 }
00888 return outputCollection;
00889 }
00890
00891 int MuonAssociatorByHits::LayerFromDetid(const DetId& detId ) const
00892 {
00893 int layerNumber=0;
00894 if (detId.det() != DetId::Tracker) return layerNumber;
00895
00896 unsigned int subdetId = static_cast<unsigned int>(detId.subdetId());
00897 if ( subdetId == StripSubdetector::TIB)
00898 {
00899 TIBDetId tibid(detId.rawId());
00900 layerNumber = tibid.layer();
00901 }
00902 else if ( subdetId == StripSubdetector::TOB )
00903 {
00904 TOBDetId tobid(detId.rawId());
00905 layerNumber = tobid.layer();
00906 }
00907 else if ( subdetId == StripSubdetector::TID)
00908 {
00909 TIDDetId tidid(detId.rawId());
00910 layerNumber = tidid.wheel();
00911 }
00912 else if ( subdetId == StripSubdetector::TEC )
00913 {
00914 TECDetId tecid(detId.rawId());
00915 layerNumber = tecid.wheel();
00916 }
00917 else if ( subdetId == PixelSubdetector::PixelBarrel )
00918 {
00919 PXBDetId pxbid(detId.rawId());
00920 layerNumber = pxbid.layer();
00921 }
00922 else if ( subdetId == PixelSubdetector::PixelEndcap )
00923 {
00924 PXFDetId pxfid(detId.rawId());
00925 layerNumber = pxfid.disk();
00926 }
00927 else edm::LogWarning("MuonAssociatorByHits")
00928 << "*** WARNING in MuonAssociatorByHits::LayerFromDetid: Unknown Tracker subdetector: subdetId = " << subdetId;
00929
00930 return layerNumber;
00931 }
00932
00933 void MuonAssociatorByHits::getMatchedIds
00934 (MapOfMatchedIds & tracker_matchedIds_valid, MapOfMatchedIds & muon_matchedIds_valid,
00935 MapOfMatchedIds & tracker_matchedIds_INVALID, MapOfMatchedIds & muon_matchedIds_INVALID,
00936 int& n_tracker_valid, int& n_dt_valid, int& n_csc_valid, int& n_rpc_valid,
00937 int& n_tracker_matched_valid, int& n_dt_matched_valid, int& n_csc_matched_valid, int& n_rpc_matched_valid,
00938 int& n_tracker_INVALID, int& n_dt_INVALID, int& n_csc_INVALID, int& n_rpc_INVALID,
00939 int& n_tracker_matched_INVALID, int& n_dt_matched_INVALID, int& n_csc_matched_INVALID, int& n_rpc_matched_INVALID,
00940 trackingRecHit_iterator begin, trackingRecHit_iterator end,
00941 TrackerHitAssociator* trackertruth,
00942 DTHitAssociator& dttruth, MuonTruth& csctruth, RPCHitAssociator& rpctruth, bool printRtS) const
00943
00944 {
00945 tracker_matchedIds_valid.clear();
00946 muon_matchedIds_valid.clear();
00947
00948 tracker_matchedIds_INVALID.clear();
00949 muon_matchedIds_INVALID.clear();
00950
00951 n_tracker_valid = 0;
00952 n_dt_valid = 0;
00953 n_csc_valid = 0;
00954 n_rpc_valid = 0;
00955
00956 n_tracker_matched_valid = 0;
00957 n_dt_matched_valid = 0;
00958 n_csc_matched_valid = 0;
00959 n_rpc_matched_valid = 0;
00960
00961 n_tracker_INVALID = 0;
00962 n_dt_INVALID = 0;
00963 n_csc_INVALID = 0;
00964 n_rpc_INVALID = 0;
00965
00966 n_tracker_matched_INVALID = 0;
00967 n_dt_matched_INVALID = 0;
00968 n_csc_matched_INVALID = 0;
00969 n_rpc_matched_INVALID = 0;
00970
00971 std::vector<SimHitIdpr> SimTrackIds;
00972
00973
00974 int iloop = 0;
00975 int iH = -1;
00976 for (trackingRecHit_iterator it = begin; it != end; it++, iloop++) {
00977 stringstream hit_index;
00978 hit_index<<iloop;
00979
00980 const TrackingRecHit * hitp = getHitPtr(it);
00981 DetId geoid = hitp->geographicalId();
00982
00983 unsigned int detid = geoid.rawId();
00984 stringstream detector_id;
00985 detector_id<<detid;
00986
00987 string hitlog = "TrackingRecHit "+hit_index.str();
00988 string wireidlog;
00989 std::vector<string> DTSimHits;
00990
00991 DetId::Detector det = geoid.det();
00992 int subdet = geoid.subdetId();
00993
00994 bool valid_Hit = hitp->isValid();
00995
00996
00997 if (det == DetId::Tracker && UseTracker) {
00998 stringstream detector_id;
00999
01000 if (subdet == PixelSubdetector::PixelBarrel) {
01001 PXBDetId pxbdetid(detid);
01002 detector_id << pxbdetid;
01003 }
01004 else if (subdet == PixelSubdetector::PixelEndcap) {
01005 PXFDetId pxfdetid(detid);
01006 detector_id << pxfdetid;
01007 }
01008 else if (subdet == StripSubdetector::TIB) {
01009 TIBDetId tibdetid(detid);
01010 detector_id << tibdetid;
01011 }
01012 else if (subdet == StripSubdetector::TOB) {
01013 TOBDetId tobdetid(detid);
01014 detector_id << tobdetid;
01015 }
01016 else if (subdet == StripSubdetector::TID) {
01017 TIDDetId tiddetid(detid);
01018 detector_id << tiddetid;
01019 }
01020 else if (subdet == StripSubdetector::TEC) {
01021 TECDetId tecdetid(detid);
01022 detector_id << tecdetid;
01023 }
01024
01025 if (valid_Hit) hitlog = hitlog+" -Tracker - detID = "+detector_id.str();
01026 else hitlog = hitlog+" *** INVALID ***"+" -Tracker - detID = "+detector_id.str();
01027
01028 iH++;
01029 SimTrackIds = trackertruth->associateHitId(*hitp);
01030
01031 if (valid_Hit) {
01032 n_tracker_valid++;
01033
01034 if(!SimTrackIds.empty()) {
01035 n_tracker_matched_valid++;
01036 tracker_matchedIds_valid[iH] = SimTrackIds;
01037 }
01038 } else {
01039 n_tracker_INVALID++;
01040
01041 if(!SimTrackIds.empty()) {
01042 n_tracker_matched_INVALID++;
01043 tracker_matchedIds_INVALID[iH] = SimTrackIds;
01044 }
01045 }
01046 }
01047
01048 else if (det == DetId::Muon && UseMuon) {
01049
01050
01051 if (subdet == MuonSubdetId::DT) {
01052 DTWireId dtdetid = DTWireId(detid);
01053 stringstream dt_detector_id;
01054 dt_detector_id << dtdetid;
01055 if (valid_Hit) hitlog = hitlog+" -Muon DT - detID = "+dt_detector_id.str();
01056 else hitlog = hitlog+" *** INVALID ***"+" -Muon DT - detID = "+dt_detector_id.str();
01057
01058 const DTRecHit1D * dtrechit = dynamic_cast<const DTRecHit1D *>(hitp);
01059
01060
01061 if (dtrechit) {
01062 iH++;
01063 SimTrackIds = dttruth.associateDTHitId(dtrechit);
01064
01065 if (valid_Hit) {
01066 n_dt_valid++;
01067
01068 if (!SimTrackIds.empty()) {
01069 n_dt_matched_valid++;
01070 muon_matchedIds_valid[iH] = SimTrackIds;
01071 }
01072 } else {
01073 n_dt_INVALID++;
01074
01075 if (!SimTrackIds.empty()) {
01076 n_dt_matched_INVALID++;
01077 muon_matchedIds_INVALID[iH] = SimTrackIds;
01078 }
01079 }
01080
01081 if (dumpDT) {
01082 DTWireId wireid = dtrechit->wireId();
01083 stringstream wid;
01084 wid<<wireid;
01085 std::vector<PSimHit> dtSimHits = dttruth.associateHit(*hitp);
01086
01087 stringstream ndthits;
01088 ndthits<<dtSimHits.size();
01089 wireidlog = "\t DTWireId :"+wid.str()+", "+ndthits.str()+" associated PSimHit :";
01090
01091 for (unsigned int j=0; j<dtSimHits.size(); j++) {
01092 stringstream index;
01093 index<<j;
01094 stringstream simhit;
01095 simhit<<dtSimHits[j];
01096 string simhitlog = "\t\t PSimHit "+index.str()+": "+simhit.str();
01097 DTSimHits.push_back(simhitlog);
01098 }
01099 }
01100 }
01101
01102
01103 else {
01104 const DTRecSegment4D * dtsegment = dynamic_cast<const DTRecSegment4D *>(hitp);
01105
01106 if (dtsegment) {
01107
01108 std::vector<const TrackingRecHit *> componentHits, phiHits, zHits;
01109 if (dtsegment->hasPhi()) {
01110 phiHits = dtsegment->phiSegment()->recHits();
01111 componentHits.insert(componentHits.end(),phiHits.begin(),phiHits.end());
01112 }
01113 if (dtsegment->hasZed()) {
01114 zHits = dtsegment->zSegment()->recHits();
01115 componentHits.insert(componentHits.end(),zHits.begin(),zHits.end());
01116 }
01117 if (printRtS) edm::LogVerbatim("MuonAssociatorByHits")
01118 <<"\n\t this TrackingRecHit is a DTRecSegment4D with "
01119 <<componentHits.size()<<" hits (phi:"<<phiHits.size()<<", z:"<<zHits.size()<<")";
01120
01121 std::vector<SimHitIdpr> i_SimTrackIds;
01122 int i_compHit = 0;
01123 for (std::vector<const TrackingRecHit *>::const_iterator ithit =componentHits.begin();
01124 ithit != componentHits.end(); ++ithit) {
01125 i_compHit++;
01126
01127 const DTRecHit1D * dtrechit1D = dynamic_cast<const DTRecHit1D *>(*ithit);
01128
01129 i_SimTrackIds.clear();
01130 if (dtrechit1D) {
01131 iH++;
01132 i_SimTrackIds = dttruth.associateDTHitId(dtrechit1D);
01133
01134 if (valid_Hit) {
01135
01136 n_dt_valid++;
01137
01138 if (!i_SimTrackIds.empty()) {
01139 n_dt_matched_valid++;
01140 muon_matchedIds_valid[iH] = i_SimTrackIds;
01141 }
01142 } else {
01143 n_dt_INVALID++;
01144
01145 if (!i_SimTrackIds.empty()) {
01146 n_dt_matched_INVALID++;
01147 muon_matchedIds_INVALID[iH] = i_SimTrackIds;
01148 }
01149 }
01150 } else if (printRtS) edm::LogWarning("MuonAssociatorByHits")
01151 <<"*** WARNING in MuonAssociatorByHits::getMatchedIds, null dynamic_cast of a DT TrackingRecHit !";
01152
01153 unsigned int i_detid = (*ithit)->geographicalId().rawId();
01154 DTWireId i_dtdetid = DTWireId(i_detid);
01155
01156 stringstream i_dt_detector_id;
01157 i_dt_detector_id << i_dtdetid;
01158
01159 stringstream i_ss;
01160 i_ss<<"\t\t hit "<<i_compHit<<" -Muon DT - detID = "<<i_dt_detector_id.str();
01161
01162 string i_hitlog = i_ss.str();
01163 i_hitlog = i_hitlog + write_matched_simtracks(i_SimTrackIds);
01164 if (printRtS) edm::LogVerbatim("MuonAssociatorByHits") << i_hitlog;
01165
01166 SimTrackIds.insert(SimTrackIds.end(),i_SimTrackIds.begin(),i_SimTrackIds.end());
01167 }
01168 }
01169
01170 else if (printRtS) edm::LogWarning("MuonAssociatorByHits")
01171 <<"*** WARNING in MuonAssociatorByHits::getMatchedIds, DT TrackingRecHit is neither DTRecHit1D nor DTRecSegment4D ! ";
01172 }
01173 }
01174
01175
01176 else if (subdet == MuonSubdetId::CSC) {
01177 CSCDetId cscdetid = CSCDetId(detid);
01178 stringstream csc_detector_id;
01179 csc_detector_id << cscdetid;
01180 if (valid_Hit) hitlog = hitlog+" -Muon CSC- detID = "+csc_detector_id.str();
01181 else hitlog = hitlog+" *** INVALID ***"+" -Muon CSC- detID = "+csc_detector_id.str();
01182
01183 const CSCRecHit2D * cscrechit = dynamic_cast<const CSCRecHit2D *>(hitp);
01184
01185
01186 if (cscrechit) {
01187 iH++;
01188 SimTrackIds = csctruth.associateCSCHitId(cscrechit);
01189
01190 if (valid_Hit) {
01191 n_csc_valid++;
01192
01193 if (!SimTrackIds.empty()) {
01194 n_csc_matched_valid++;
01195 muon_matchedIds_valid[iH] = SimTrackIds;
01196 }
01197 } else {
01198 n_csc_INVALID++;
01199
01200 if (!SimTrackIds.empty()) {
01201 n_csc_matched_INVALID++;
01202 muon_matchedIds_INVALID[iH] = SimTrackIds;
01203 }
01204 }
01205 }
01206
01207
01208 else {
01209 const CSCSegment * cscsegment = dynamic_cast<const CSCSegment *>(hitp);
01210
01211 if (cscsegment) {
01212
01213 std::vector<const TrackingRecHit *> componentHits = cscsegment->recHits();
01214 if (printRtS) edm::LogVerbatim("MuonAssociatorByHits")
01215 <<"\n\t this TrackingRecHit is a CSCSegment with "<<componentHits.size()<<" hits";
01216
01217 std::vector<SimHitIdpr> i_SimTrackIds;
01218 int i_compHit = 0;
01219 for (std::vector<const TrackingRecHit *>::const_iterator ithit =componentHits.begin();
01220 ithit != componentHits.end(); ++ithit) {
01221 i_compHit++;
01222
01223 const CSCRecHit2D * cscrechit2D = dynamic_cast<const CSCRecHit2D *>(*ithit);
01224
01225 i_SimTrackIds.clear();
01226 if (cscrechit2D) {
01227 iH++;
01228 i_SimTrackIds = csctruth.associateCSCHitId(cscrechit2D);
01229
01230 if (valid_Hit) {
01231
01232 n_csc_valid++;
01233
01234 if (!i_SimTrackIds.empty()) {
01235 n_csc_matched_valid++;
01236 muon_matchedIds_valid[iH] = i_SimTrackIds;
01237 }
01238 } else {
01239 n_csc_INVALID++;
01240
01241 if (!i_SimTrackIds.empty()) {
01242 n_csc_matched_INVALID++;
01243 muon_matchedIds_INVALID[iH] = i_SimTrackIds;
01244 }
01245 }
01246 } else if (printRtS) edm::LogWarning("MuonAssociatorByHits")
01247 <<"*** WARNING in MuonAssociatorByHits::getMatchedIds, null dynamic_cast of a CSC TrackingRecHit !";
01248
01249 unsigned int i_detid = (*ithit)->geographicalId().rawId();
01250 CSCDetId i_cscdetid = CSCDetId(i_detid);
01251
01252 stringstream i_csc_detector_id;
01253 i_csc_detector_id << i_cscdetid;
01254
01255 stringstream i_ss;
01256 i_ss<<"\t\t hit "<<i_compHit<<" -Muon CSC- detID = "<<i_csc_detector_id.str();
01257
01258 string i_hitlog = i_ss.str();
01259 i_hitlog = i_hitlog + write_matched_simtracks(i_SimTrackIds);
01260 if (printRtS) edm::LogVerbatim("MuonAssociatorByHits") << i_hitlog;
01261
01262 SimTrackIds.insert(SimTrackIds.end(),i_SimTrackIds.begin(),i_SimTrackIds.end());
01263 }
01264 }
01265
01266 else if (printRtS) edm::LogWarning("MuonAssociatorByHits")
01267 <<"*** WARNING in MuonAssociatorByHits::getMatchedIds, CSC TrackingRecHit is neither CSCRecHit2D nor CSCSegment ! ";
01268 }
01269 }
01270
01271
01272 else if (subdet == MuonSubdetId::RPC) {
01273 RPCDetId rpcdetid = RPCDetId(detid);
01274 stringstream rpc_detector_id;
01275 rpc_detector_id << rpcdetid;
01276 if (valid_Hit) hitlog = hitlog+" -Muon RPC- detID = "+rpc_detector_id.str();
01277 else hitlog = hitlog+" *** INVALID ***"+" -Muon RPC- detID = "+rpc_detector_id.str();
01278
01279 iH++;
01280 SimTrackIds = rpctruth.associateRecHit(*hitp);
01281
01282 if (valid_Hit) {
01283 n_rpc_valid++;
01284
01285 if (!SimTrackIds.empty()) {
01286 n_rpc_matched_valid++;
01287 muon_matchedIds_valid[iH] = SimTrackIds;
01288 }
01289 } else {
01290 n_rpc_INVALID++;
01291
01292 if (!SimTrackIds.empty()) {
01293 n_rpc_matched_INVALID++;
01294 muon_matchedIds_INVALID[iH] = SimTrackIds;
01295 }
01296 }
01297
01298 } else if (printRtS) edm::LogVerbatim("MuonAssociatorByHits")
01299 <<"TrackingRecHit "<<iloop<<" *** WARNING *** Unexpected Hit from Detector = "<<det;
01300 }
01301 else continue;
01302
01303 hitlog = hitlog + write_matched_simtracks(SimTrackIds);
01304
01305 if (printRtS) edm::LogVerbatim("MuonAssociatorByHits") << hitlog;
01306 if (printRtS && dumpDT && det==DetId::Muon && subdet==MuonSubdetId::DT) {
01307 edm::LogVerbatim("MuonAssociatorByHits") <<wireidlog;
01308 for (unsigned int j=0; j<DTSimHits.size(); j++) {
01309 edm::LogVerbatim("MuonAssociatorByHits") <<DTSimHits[j];
01310 }
01311 }
01312
01313 }
01314 }
01315
01316 int MuonAssociatorByHits::getShared(MapOfMatchedIds & matchedIds, TrackingParticleCollection::const_iterator trpart) const {
01317 int nshared = 0;
01318
01319
01320 for (MapOfMatchedIds::const_iterator iRecH=matchedIds.begin(); iRecH!=matchedIds.end(); ++iRecH) {
01321
01322
01323 std::vector<SimHitIdpr> SimTrackIds = (*iRecH).second;
01324
01325 bool found = false;
01326
01327 for (std::vector<SimHitIdpr>::const_iterator iSimH=SimTrackIds.begin(); iSimH!=SimTrackIds.end(); ++iSimH) {
01328 uint32_t simtrackId = iSimH->first;
01329 EncodedEventId evtId = iSimH->second;
01330
01331
01332 for (TrackingParticle::g4t_iterator simtrack = trpart->g4Track_begin(); simtrack != trpart->g4Track_end(); ++simtrack) {
01333 if (simtrack->trackId() == simtrackId && simtrack->eventId() == evtId) {
01334 found = true;
01335 break;
01336 }
01337 }
01338
01339 if (found) {
01340 nshared++;
01341 break;
01342 }
01343 }
01344 }
01345
01346 return nshared;
01347 }
01348
01349 std::string MuonAssociatorByHits::write_matched_simtracks(const std::vector<SimHitIdpr>& SimTrackIds) const {
01350
01351 string hitlog;
01352
01353 if (!SimTrackIds.empty()) {
01354 hitlog = " matched to SimTrack";
01355
01356 for(size_t j=0; j<SimTrackIds.size(); j++){
01357 stringstream trackid;
01358 trackid<<SimTrackIds[j].first;
01359
01360 stringstream evtid;
01361 evtid<<SimTrackIds[j].second.event();
01362
01363 stringstream bunchxid;
01364 bunchxid<<SimTrackIds[j].second.bunchCrossing();
01365
01366 hitlog = hitlog+" Id:"+trackid.str()+"/Evt:("+evtid.str()+","+bunchxid.str()+") ";
01367 }
01368 } else hitlog = " *** UNMATCHED ***";
01369
01370 return hitlog;
01371 }
01372
01373 void MuonAssociatorByHits::associateMuons(MuonToSimCollection & recToSim, SimToMuonCollection & simToRec,
01374 const edm::Handle<edm::View<reco::Muon> > &tCH, MuonTrackType type,
01375 const edm::Handle<TrackingParticleCollection>&tPCH,
01376 const edm::Event * event, const edm::EventSetup * setup) const {
01377
01378 edm::RefVector<TrackingParticleCollection> tpc(tPCH.id());
01379 for (unsigned int j=0; j<tPCH->size();j++)
01380 tpc.push_back(edm::Ref<TrackingParticleCollection>(tPCH,j));
01381
01382 associateMuons(recToSim, simToRec, tCH->refVector(),type,tpc,event,setup);
01383 }
01384
01385 void MuonAssociatorByHits::associateMuons(MuonToSimCollection & recToSim, SimToMuonCollection & simToRec,
01386 const edm::RefToBaseVector<reco::Muon> & muons, MuonTrackType trackType,
01387 const edm::RefVector<TrackingParticleCollection>& tPC,
01388 const edm::Event * event, const edm::EventSetup * setup) const {
01389
01391 MuonAssociatorByHits::TrackHitsCollection muonHitRefs;
01392 edm::OwnVector<TrackingRecHit> allTMRecHits;
01393 TrackingRecHitRefVector hitRefVector;
01394 switch (trackType) {
01395 case InnerTk:
01396 for (edm::RefToBaseVector<reco::Muon>::const_iterator it = muons.begin(), ed = muons.end(); it != ed; ++it) {
01397 edm::RefToBase<reco::Muon> mur = *it;
01398 if (mur->track().isNonnull()) {
01399 muonHitRefs.push_back(std::make_pair(mur->track()->recHitsBegin(), mur->track()->recHitsEnd()));
01400 } else {
01401 muonHitRefs.push_back(std::make_pair(hitRefVector.begin(), hitRefVector.end()));
01402 }
01403 }
01404 break;
01405 case OuterTk:
01406 for (edm::RefToBaseVector<reco::Muon>::const_iterator it = muons.begin(), ed = muons.end(); it != ed; ++it) {
01407 edm::RefToBase<reco::Muon> mur = *it;
01408 if (mur->outerTrack().isNonnull()) {
01409 muonHitRefs.push_back(std::make_pair(mur->outerTrack()->recHitsBegin(), mur->outerTrack()->recHitsEnd()));
01410 } else {
01411 muonHitRefs.push_back(std::make_pair(hitRefVector.begin(), hitRefVector.end()));
01412 }
01413 }
01414 break;
01415 case GlobalTk:
01416 for (edm::RefToBaseVector<reco::Muon>::const_iterator it = muons.begin(), ed = muons.end(); it != ed; ++it) {
01417 edm::RefToBase<reco::Muon> mur = *it;
01418 if (mur->globalTrack().isNonnull()) {
01419 muonHitRefs.push_back(std::make_pair(mur->globalTrack()->recHitsBegin(), mur->globalTrack()->recHitsEnd()));
01420 } else {
01421 muonHitRefs.push_back(std::make_pair(hitRefVector.begin(), hitRefVector.end()));
01422 }
01423 }
01424 break;
01425 case Segments: {
01426 TrackerMuonHitExtractor hitExtractor(conf_);
01427 hitExtractor.init(*event, *setup);
01428
01429 std::vector<std::pair<size_t, size_t> > muonHitIndices;
01430 for (edm::RefToBaseVector<reco::Muon>::const_iterator it = muons.begin(), ed = muons.end(); it != ed; ++it) {
01431 edm::RefToBase<reco::Muon> mur = *it;
01432 std::pair<size_t, size_t> indices(allTMRecHits.size(), allTMRecHits.size());
01433 if (mur->isTrackerMuon()) {
01434 std::vector<const TrackingRecHit *> hits = hitExtractor.getMuonHits(*mur);
01435 for (std::vector<const TrackingRecHit *>::const_iterator ith = hits.begin(), edh = hits.end(); ith != edh; ++ith) {
01436 allTMRecHits.push_back(**ith);
01437 }
01438 indices.second += hits.size();
01439 }
01440 muonHitIndices.push_back(indices);
01441 }
01442
01443 for (size_t i = 0, n = allTMRecHits.size(); i < n; ++i) {
01444 hitRefVector.push_back(TrackingRecHitRef(& allTMRecHits, i));
01445 }
01446
01447 typedef std::pair<size_t, size_t> index_pair;
01448 trackingRecHit_iterator hitRefBegin = hitRefVector.begin();
01449 for (std::vector<std::pair<size_t, size_t> >::const_iterator idxs = muonHitIndices.begin(), idxend = muonHitIndices.end(); idxs != idxend; ++idxs) {
01450 muonHitRefs.push_back(std::make_pair(hitRefBegin+idxs->first,
01451 hitRefBegin+idxs->second));
01452 }
01453
01454 }
01455 break;
01456 }
01457
01459 MuonAssociatorByHits::IndexAssociation recSimColl = associateRecoToSimIndices(muonHitRefs,tPC,event,setup);
01460 for (MuonAssociatorByHits::IndexAssociation::const_iterator it = recSimColl.begin(), ed = recSimColl.end(); it != ed; ++it) {
01461 edm::RefToBase<reco::Muon> rec = muons[it->first];
01462 const std::vector<MuonAssociatorByHits::IndexMatch> & idxAss = it->second;
01463 std::vector<std::pair<TrackingParticleRef, double> > & tpAss = recToSim[rec];
01464 for (std::vector<MuonAssociatorByHits::IndexMatch>::const_iterator ita = idxAss.begin(), eda = idxAss.end(); ita != eda; ++ita) {
01465 tpAss.push_back(std::make_pair(tPC[ita->idx], ita->quality));
01466 }
01467 }
01468 MuonAssociatorByHits::IndexAssociation simRecColl = associateSimToRecoIndices(muonHitRefs,tPC,event,setup);
01469 for (MuonAssociatorByHits::IndexAssociation::const_iterator it = simRecColl.begin(), ed = simRecColl.end(); it != ed; ++it) {
01470 TrackingParticleRef sim = tPC[it->first];
01471 const std::vector<MuonAssociatorByHits::IndexMatch> & idxAss = it->second;
01472 std::vector<std::pair<edm::RefToBase<reco::Muon>, double> > & recAss = simToRec[sim];
01473 for (std::vector<MuonAssociatorByHits::IndexMatch>::const_iterator ita = idxAss.begin(), eda = idxAss.end(); ita != eda; ++ita) {
01474 recAss.push_back(std::make_pair(muons[ita->idx], ita->quality));
01475 }
01476 }
01477
01478 }
1.7.3