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