Reference Manual

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