CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/SimMuon/MCTruth/src/MuonAssociatorByHits.cc

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_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     // all used hits (valid+INVALID), defined by UseTracker, UseMuon
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     // selected hits are set initially to valid hits
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     // matched hits are a subsample of the selected hits
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       // selected muon hits = INVALID when (useZeroHitMuons == True) and track has no valid muon hits
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         // global_quality used to order the matching TPs
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           //if a track has just 3 hits in the Tracker we require that all 3 hits are shared
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         // (matchOk) has to account for different track types (tracker-only, standalone muons, global muons)
00396         bool matchOk = trackerOk || muonOk;
00397 
00398         // only for global muons: match both tracker and muon stub unless (acceptOneStubMatchings==true)
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           // print something only if this TrackingParticle shares some hits with the current reco::Track
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       }    //  loop over TrackingParticle
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     }    //  if(n_matching_simhits != 0)
00441     
00442   }    // loop over reco::Track
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(); // perhaps not even necessary
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   // Tracker hit association  
00512   TrackerHitAssociator * trackertruth = new TrackerHitAssociator(*e, conf_);
00513   // CSC hit association
00514   MuonTruth csctruth(*e,*setup,conf_);
00515   // DT hit association
00516   printRtS = false;
00517   DTHitAssociator dttruth(*e,*setup,conf_,printRtS);  
00518   // RPC hit association
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     // all hits = valid +INVALID
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     // all used hits (valid+INVALID), defined by UseTracker, UseMuon
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      // selected hits are set initially to valid hits
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     // matched hits are a subsample of the selected hits
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       // selected muon hits = INVALID when (useZeroHitMuons == True) and track has no valid muon hits
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         // shared hits are counted over the selected ones
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; // if this TP shares no hits with the current reco::Track loop over 
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               //no grouped, no splitting
00701               if (!UseGrouped && !UseSplitting)
00702                 if (LayerFromDetid(dId)==LayerFromDetid(dIdOK) &&
00703                     dId.subdetId()==dIdOK.subdetId()) newhit = false;
00704               //no grouped, splitting
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               //grouped, no splitting
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               //grouped, splitting
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             // discard BAD CSC chambers (ME4/2) from hit counting
00729             if (dId.subdetId() == MuonSubdetId::CSC) {
00730               if (csctruth.cscBadChambers->isInBadChamber(CSCDetId(dId))) {
00731                 // edm::LogVerbatim("MuonAssociatorByHits")<<"This PSimHit is in a BAD CSC chamber, CSCDetId = "<<CSCDetId(dId);
00732                 n_muon_simhits--;
00733               }
00734             }
00735             
00736           }
00737         }
00738 
00739         n_tracker_recounted_simhits = tphits.size();
00740         // Handle the case of TrackingParticles that don't have PSimHits inside, e.g. because they were made on RECOSIM only.
00741         if (trpart->trackPSimHit().empty()) {
00742             // FIXME this can be made better, counting the digiSimLinks associated to this TP, but perhaps it's not worth it
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         // global_quality used to order the matching tracks
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         // global purity
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           //if a track has just 3 hits in the Tracker we require that all 3 hits are shared
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         // (matchOk) has to account for different track types (tracker-only, standalone muons, global muons)
00809         bool matchOk = trackerOk || muonOk;
00810 
00811         // only for global muons: match both tracker and muon stub unless (acceptOneStubMatchings==true)
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           // print something only if this TrackingParticle shares some hits with the current reco::Track
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       }  // loop over TrackingParticle's
00870     }   // if(n_matching_simhits != 0)
00871   }    // loop over reco Tracks
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   // main loop on TrackingRecHits
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     // Si-Tracker Hits
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     // Muon detector Hits    
01048     else if (det == DetId::Muon && UseMuon) {
01049 
01050       // DT Hits      
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         // single DT hits
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           }  // if (dumpDT)
01100         }
01101 
01102         // DT segments  
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                   // validity check is on the segment, but hits are counted one-by-one
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           }  // if (dtsegment)
01169 
01170           else if (printRtS) edm::LogWarning("MuonAssociatorByHits")
01171             <<"*** WARNING in MuonAssociatorByHits::getMatchedIds, DT TrackingRecHit is neither DTRecHit1D nor DTRecSegment4D ! ";          
01172         }
01173       }
01174       
01175       // CSC Hits
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         // single CSC hits
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         // CSC segments
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                   // validity check is on the segment, but hits are counted one-by-one
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           }  // if (cscsegment)
01265 
01266           else if (printRtS) edm::LogWarning("MuonAssociatorByHits")
01267             <<"*** WARNING in MuonAssociatorByHits::getMatchedIds, CSC TrackingRecHit is neither CSCRecHit2D nor CSCSegment ! ";
01268         }
01269       }
01270       
01271       // RPC Hits
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   } //trackingRecHit loop
01314 }
01315 
01316 int MuonAssociatorByHits::getShared(MapOfMatchedIds & matchedIds, TrackingParticleCollection::const_iterator trpart) const {
01317   int nshared = 0;
01318 
01319   // map is indexed over the rechits of the reco::Track (no double-countings allowed)
01320   for (MapOfMatchedIds::const_iterator iRecH=matchedIds.begin(); iRecH!=matchedIds.end(); ++iRecH) {
01321 
01322     // vector of associated simhits associated to the current rechit
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       // look for shared hits with the given TrackingParticle (looping over component SimTracks)
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;  // this I will fill in only for tracker muon hits from segments
01393     TrackingRecHitRefVector hitRefVector;              // same as above, plus used to get null iterators for muons without a track
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                 // puts hits in the vector, and record indices
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                 // puts hits in the ref-vector
01443                 for (size_t i = 0, n = allTMRecHits.size(); i < n; ++i) {
01444                     hitRefVector.push_back(TrackingRecHitRef(& allTMRecHits, i)); 
01445                 }
01446                 // convert indices into pairs of iterators to references
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 }
Generated for CMSSW Reference Manual by doxygen 1.7.3