Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-10 02:59:12 

0001 #include "DataFormats/CSCRecHit/interface/CSCSegment.h"
0002 #include "DataFormats/DTRecHit/interface/DTRecSegment4D.h"
0003 #include "DataFormats/DetId/interface/DetId.h"
0004 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0005 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0006 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0007 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0010 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0011 #include "SimMuon/MCTruth/interface/MuonAssociatorByHitsHelper.h"
0012 #include <sstream>
0013 
0014 using namespace reco;
0015 using namespace std;
0016 
0017 MuonAssociatorByHitsHelper::MuonAssociatorByHitsHelper(const edm::ParameterSet &conf)
0018     : includeZeroHitMuons(conf.getParameter<bool>("includeZeroHitMuons")),
0019       acceptOneStubMatchings(conf.getParameter<bool>("acceptOneStubMatchings")),
0020       rejectBadGlobal(conf.getParameter<bool>("rejectBadGlobal")),
0021       UseTracker(conf.getParameter<bool>("UseTracker")),
0022       UseMuon(conf.getParameter<bool>("UseMuon")),
0023       AbsoluteNumberOfHits_track(conf.getParameter<bool>("AbsoluteNumberOfHits_track")),
0024       NHitCut_track(conf.getParameter<unsigned int>("NHitCut_track")),
0025       EfficiencyCut_track(conf.getParameter<double>("EfficiencyCut_track")),
0026       PurityCut_track(conf.getParameter<double>("PurityCut_track")),
0027       AbsoluteNumberOfHits_muon(conf.getParameter<bool>("AbsoluteNumberOfHits_muon")),
0028       NHitCut_muon(conf.getParameter<unsigned int>("NHitCut_muon")),
0029       EfficiencyCut_muon(conf.getParameter<double>("EfficiencyCut_muon")),
0030       PurityCut_muon(conf.getParameter<double>("PurityCut_muon")),
0031       UsePixels(conf.getParameter<bool>("UsePixels")),
0032       UseGrouped(conf.getParameter<bool>("UseGrouped")),
0033       UseSplitting(conf.getParameter<bool>("UseSplitting")),
0034       ThreeHitTracksAreSpecial(conf.getParameter<bool>("ThreeHitTracksAreSpecial")),
0035       dumpDT(conf.getParameter<bool>("dumpDT")) {
0036   edm::LogVerbatim("MuonAssociatorByHitsHelper") << "\n constructing  MuonAssociatorByHitsHelper" << conf.dump();
0037 
0038   // up to the user in the other cases - print a message
0039   if (UseTracker)
0040     edm::LogVerbatim("MuonAssociatorByHitsHelper") << "\n UseTracker = TRUE  : Tracker SimHits and RecHits WILL be "
0041                                                       "counted";
0042   else
0043     edm::LogVerbatim("MuonAssociatorByHitsHelper") << "\n UseTracker = FALSE : Tracker SimHits and RecHits WILL NOT be "
0044                                                       "counted";
0045 
0046   // up to the user in the other cases - print a message
0047   if (UseMuon)
0048     edm::LogVerbatim("MuonAssociatorByHitsHelper") << " UseMuon = TRUE  : Muon SimHits and RecHits WILL be counted";
0049   else
0050     edm::LogVerbatim("MuonAssociatorByHitsHelper")
0051         << " UseMuon = FALSE : Muon SimHits and RecHits WILL NOT be counted" << endl;
0052 
0053   // check consistency of the configuration when allowing zero-hit muon matching
0054   // (counting invalid hits)
0055   if (includeZeroHitMuons) {
0056     edm::LogVerbatim("MuonAssociatorByHitsHelper") << "\n includeZeroHitMuons = TRUE"
0057                                                    << "\n ==> (re)set NHitCut_muon = 0, PurityCut_muon = 0, "
0058                                                       "EfficiencyCut_muon = 0"
0059                                                    << endl;
0060     NHitCut_muon = 0;
0061     PurityCut_muon = 0.;
0062     EfficiencyCut_muon = 0.;
0063   }
0064 }
0065 
0066 MuonAssociatorByHitsHelper::IndexAssociation MuonAssociatorByHitsHelper::associateRecoToSimIndices(
0067     const TrackHitsCollection &tC,
0068     const edm::RefVector<TrackingParticleCollection> &TPCollectionH,
0069     const Resources &resources) const {
0070   auto tTopo = resources.tTopo_;
0071   auto trackertruth = resources.trackerHitAssoc_;
0072   auto const &csctruth = *resources.cscHitAssoc_;
0073   auto const &dttruth = *resources.dtHitAssoc_;
0074   auto const &rpctruth = *resources.rpcHitAssoc_;
0075   auto const &gemtruth = *resources.gemHitAssoc_;
0076 
0077   int tracker_nshared = 0;
0078   int muon_nshared = 0;
0079   int global_nshared = 0;
0080 
0081   double tracker_quality = 0;
0082   double tracker_quality_cut;
0083   if (AbsoluteNumberOfHits_track)
0084     tracker_quality_cut = static_cast<double>(NHitCut_track);
0085   else
0086     tracker_quality_cut = PurityCut_track;
0087 
0088   double muon_quality = 0;
0089   double muon_quality_cut;
0090   if (AbsoluteNumberOfHits_muon)
0091     muon_quality_cut = static_cast<double>(NHitCut_muon);
0092   else
0093     muon_quality_cut = PurityCut_muon;
0094 
0095   double global_quality = 0;
0096 
0097   MapOfMatchedIds tracker_matchedIds_valid, muon_matchedIds_valid;
0098   MapOfMatchedIds tracker_matchedIds_INVALID, muon_matchedIds_INVALID;
0099 
0100   IndexAssociation outputCollection;
0101 
0102   TrackingParticleCollection tPC;
0103   tPC.reserve(TPCollectionH.size());
0104   for (auto const &ref : TPCollectionH) {
0105     tPC.push_back(*ref);
0106   }
0107 
0108   if (resources.diagnostics_) {
0109     resources.diagnostics_(tC, tPC);
0110   }
0111 
0112   int tindex = 0;
0113   for (TrackHitsCollection::const_iterator track = tC.begin(); track != tC.end(); track++, tindex++) {
0114     edm::LogVerbatim("MuonAssociatorByHitsHelper")
0115         << "\n"
0116         << "reco::Track " << tindex << ", number of RecHits = " << (track->second - track->first) << "\n";
0117     tracker_matchedIds_valid.clear();
0118     muon_matchedIds_valid.clear();
0119 
0120     tracker_matchedIds_INVALID.clear();
0121     muon_matchedIds_INVALID.clear();
0122 
0123     bool this_track_matched = false;
0124     int n_matching_simhits = 0;
0125 
0126     // all hits = valid +INVALID
0127     int n_all = 0;
0128     int n_tracker_all = 0;
0129     int n_dt_all = 0;
0130     int n_csc_all = 0;
0131     int n_rpc_all = 0;
0132     int n_gem_all = 0;
0133 
0134     int n_valid = 0;
0135     int n_tracker_valid = 0;
0136     int n_muon_valid = 0;
0137     int n_dt_valid = 0;
0138     int n_csc_valid = 0;
0139     int n_rpc_valid = 0;
0140     int n_gem_valid = 0;
0141 
0142     int n_tracker_matched_valid = 0;
0143     int n_muon_matched_valid = 0;
0144     int n_dt_matched_valid = 0;
0145     int n_csc_matched_valid = 0;
0146     int n_rpc_matched_valid = 0;
0147     int n_gem_matched_valid = 0;
0148 
0149     int n_INVALID = 0;
0150     int n_tracker_INVALID = 0;
0151     int n_muon_INVALID = 0;
0152     int n_dt_INVALID = 0;
0153     int n_csc_INVALID = 0;
0154     int n_rpc_INVALID = 0;
0155     int n_gem_INVALID = 0;
0156 
0157     int n_tracker_matched_INVALID = 0;
0158     int n_muon_matched_INVALID = 0;
0159     int n_dt_matched_INVALID = 0;
0160     int n_csc_matched_INVALID = 0;
0161     int n_rpc_matched_INVALID = 0;
0162     int n_gem_matched_INVALID = 0;
0163 
0164     bool printRtS = true;
0165     getMatchedIds(tracker_matchedIds_valid,
0166                   muon_matchedIds_valid,
0167                   tracker_matchedIds_INVALID,
0168                   muon_matchedIds_INVALID,
0169                   n_tracker_valid,
0170                   n_dt_valid,
0171                   n_csc_valid,
0172                   n_rpc_valid,
0173                   n_gem_valid,
0174                   n_tracker_matched_valid,
0175                   n_dt_matched_valid,
0176                   n_csc_matched_valid,
0177                   n_rpc_matched_valid,
0178                   n_gem_matched_valid,
0179                   n_tracker_INVALID,
0180                   n_dt_INVALID,
0181                   n_csc_INVALID,
0182                   n_rpc_INVALID,
0183                   n_gem_INVALID,
0184                   n_tracker_matched_INVALID,
0185                   n_dt_matched_INVALID,
0186                   n_csc_matched_INVALID,
0187                   n_rpc_matched_INVALID,
0188                   n_gem_matched_INVALID,
0189                   track->first,
0190                   track->second,
0191                   trackertruth,
0192                   dttruth,
0193                   csctruth,
0194                   rpctruth,
0195                   gemtruth,
0196                   printRtS,
0197                   tTopo);
0198 
0199     n_matching_simhits = tracker_matchedIds_valid.size() + muon_matchedIds_valid.size() +
0200                          tracker_matchedIds_INVALID.size() + muon_matchedIds_INVALID.size();
0201 
0202     n_muon_valid = n_dt_valid + n_csc_valid + n_rpc_valid + n_gem_valid;
0203     n_valid = n_tracker_valid + n_muon_valid;
0204     n_muon_INVALID = n_dt_INVALID + n_csc_INVALID + n_rpc_INVALID + n_gem_INVALID;
0205     n_INVALID = n_tracker_INVALID + n_muon_INVALID;
0206 
0207     // all used hits (valid+INVALID), defined by UseTracker, UseMuon
0208     n_tracker_all = n_tracker_valid + n_tracker_INVALID;
0209     n_dt_all = n_dt_valid + n_dt_INVALID;
0210     n_csc_all = n_csc_valid + n_csc_INVALID;
0211     n_rpc_all = n_rpc_valid + n_rpc_INVALID;
0212     n_gem_all = n_gem_valid + n_gem_INVALID;
0213     n_all = n_valid + n_INVALID;
0214 
0215     n_muon_matched_valid = n_dt_matched_valid + n_csc_matched_valid + n_rpc_matched_valid + n_gem_matched_valid;
0216     n_muon_matched_INVALID =
0217         n_dt_matched_INVALID + n_csc_matched_INVALID + n_rpc_matched_INVALID + n_gem_matched_INVALID;
0218 
0219     // selected hits are set initially to valid hits
0220     int n_tracker_selected_hits = n_tracker_valid;
0221     int n_muon_selected_hits = n_muon_valid;
0222     int n_dt_selected_hits = n_dt_valid;
0223     int n_csc_selected_hits = n_csc_valid;
0224     int n_rpc_selected_hits = n_rpc_valid;
0225     int n_gem_selected_hits = n_gem_valid;
0226 
0227     // matched hits are a subsample of the selected hits
0228     int n_tracker_matched = n_tracker_matched_valid;
0229     int n_muon_matched = n_muon_matched_valid;
0230     int n_dt_matched = n_dt_matched_valid;
0231     int n_csc_matched = n_csc_matched_valid;
0232     int n_rpc_matched = n_rpc_matched_valid;
0233     int n_gem_matched = n_gem_matched_valid;
0234 
0235     std::string InvMuonHits, ZeroHitMuon;
0236 
0237     if (includeZeroHitMuons && n_muon_valid == 0 && n_muon_INVALID != 0) {
0238       // selected muon hits = INVALID when (useZeroHitMuons == True) and track
0239       // has no valid muon hits
0240 
0241       InvMuonHits = " ***INVALID MUON HITS***";
0242       ZeroHitMuon = " ***ZERO-HIT MUON***";
0243 
0244       n_muon_selected_hits = n_muon_INVALID;
0245       n_dt_selected_hits = n_dt_INVALID;
0246       n_csc_selected_hits = n_csc_INVALID;
0247       n_rpc_selected_hits = n_rpc_INVALID;
0248       n_gem_selected_hits = n_gem_INVALID;
0249 
0250       n_muon_matched = n_muon_matched_INVALID;
0251       n_dt_matched = n_dt_matched_INVALID;
0252       n_csc_matched = n_csc_matched_INVALID;
0253       n_rpc_matched = n_rpc_matched_INVALID;
0254       n_gem_matched = n_gem_matched_INVALID;
0255     }
0256 
0257     int n_selected_hits = n_tracker_selected_hits + n_muon_selected_hits;
0258     int n_matched = n_tracker_matched + n_muon_matched;
0259 
0260     edm::LogVerbatim("MuonAssociatorByHitsHelper")
0261         << "\n"
0262         << "# TrackingRecHits: " << (track->second - track->first) << "\n"
0263         << "# used RecHits     = " << n_all << " (" << n_tracker_all << "/" << n_dt_all << "/" << n_csc_all << "/"
0264         << n_rpc_all << "/" << n_gem_all << " in Tracker/DT/CSC/RPC/GEM)"
0265         << ", obtained from " << n_matching_simhits << " SimHits"
0266         << "\n"
0267         << "# selected RecHits = " << n_selected_hits << " (" << n_tracker_selected_hits << "/" << n_dt_selected_hits
0268         << "/" << n_csc_selected_hits << "/" << n_rpc_selected_hits << "/" << n_gem_selected_hits
0269         << " in Tracker/DT/CSC/RPC/GEM)" << InvMuonHits << "\n"
0270         << "# matched RecHits  = " << n_matched << " (" << n_tracker_matched << "/" << n_dt_matched << "/"
0271         << n_csc_matched << "/" << n_rpc_matched << "/" << n_gem_matched << " in Tracker/DT/CSC/RPC/GEM)";
0272 
0273     if (n_all > 0 && n_matching_simhits == 0)
0274       edm::LogVerbatim("MuonAssociatorByHitsHelper")
0275           << "*** WARNING in MuonAssociatorByHitsHelper::associateRecoToSim: "
0276              "no matching PSimHit found for this reco::Track !";
0277 
0278     if (n_matching_simhits != 0) {
0279       edm::LogVerbatim("MuonAssociatorByHitsHelper")
0280           << "\n"
0281           << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
0282              "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
0283           << "\n"
0284           << "reco::Track " << tindex << ZeroHitMuon << "\n\t"
0285           << "made of " << n_selected_hits << " selected RecHits (tracker:" << n_tracker_selected_hits
0286           << "/muons:" << n_muon_selected_hits << ")";
0287 
0288       int tpindex = 0;
0289       for (TrackingParticleCollection::const_iterator trpart = tPC.begin(); trpart != tPC.end(); ++trpart, ++tpindex) {
0290         tracker_nshared = getShared(tracker_matchedIds_valid, trpart);
0291         muon_nshared = getShared(muon_matchedIds_valid, trpart);
0292 
0293         if (includeZeroHitMuons && n_muon_valid == 0 && n_muon_INVALID != 0)
0294           muon_nshared = getShared(muon_matchedIds_INVALID, trpart);
0295 
0296         global_nshared = tracker_nshared + muon_nshared;
0297 
0298         if (AbsoluteNumberOfHits_track)
0299           tracker_quality = static_cast<double>(tracker_nshared);
0300         else if (n_tracker_selected_hits != 0)
0301           tracker_quality = (static_cast<double>(tracker_nshared) / static_cast<double>(n_tracker_selected_hits));
0302         else
0303           tracker_quality = 0;
0304 
0305         if (AbsoluteNumberOfHits_muon)
0306           muon_quality = static_cast<double>(muon_nshared);
0307         else if (n_muon_selected_hits != 0)
0308           muon_quality = (static_cast<double>(muon_nshared) / static_cast<double>(n_muon_selected_hits));
0309         else
0310           muon_quality = 0;
0311 
0312         // global_quality used to order the matching TPs
0313         if (n_selected_hits != 0) {
0314           if (AbsoluteNumberOfHits_muon && AbsoluteNumberOfHits_track)
0315             global_quality = global_nshared;
0316           else
0317             global_quality = (static_cast<double>(global_nshared) / static_cast<double>(n_selected_hits));
0318         } else
0319           global_quality = 0;
0320 
0321         bool trackerOk = false;
0322         if (n_tracker_selected_hits != 0) {
0323           if (tracker_quality > tracker_quality_cut)
0324             trackerOk = true;
0325           // if a track has just 3 hits in the Tracker we require that all 3
0326           // hits are shared
0327           if (ThreeHitTracksAreSpecial && n_tracker_selected_hits == 3 && tracker_nshared < 3)
0328             trackerOk = false;
0329         }
0330 
0331         bool muonOk = false;
0332         if (n_muon_selected_hits != 0) {
0333           if (muon_quality > muon_quality_cut)
0334             muonOk = true;
0335         }
0336 
0337         // (matchOk) has to account for different track types (tracker-only,
0338         // standalone muons, global muons)
0339         bool matchOk = trackerOk || muonOk;
0340 
0341         // only for global tracks: match both tracker and muon stub (if acceptOneStubMatchings==false)
0342         // depending on the muon selection reject tracks with only one stub (if rejectBadGlobal==true)
0343         //
0344         if (UseTracker && UseMuon && !acceptOneStubMatchings &&
0345             ((n_tracker_selected_hits != 0 && n_muon_selected_hits != 0) || rejectBadGlobal))
0346           matchOk = trackerOk && muonOk;
0347 
0348         if (matchOk) {
0349           outputCollection[tindex].push_back(IndexMatch(tpindex, global_quality));
0350           this_track_matched = true;
0351 
0352           edm::LogVerbatim("MuonAssociatorByHitsHelper")
0353               << "\n\t"
0354               << " **MATCHED** with quality = " << global_quality << " (tracker: " << tracker_quality
0355               << " / muon: " << muon_quality << ")"
0356               << "\n\t"
0357               << "   N shared hits = " << global_nshared << " (tracker: " << tracker_nshared
0358               << " / muon: " << muon_nshared << ")"
0359               << "\n"
0360               << "   to TrackingParticle:  q = " << (*trpart).charge() << ", p = " << (*trpart).p()
0361               << ", pT = " << (*trpart).pt() << ", eta = " << (*trpart).eta() << ", phi = " << (*trpart).phi() << "\n\t"
0362               << " pdg code = " << (*trpart).pdgId() << ", made of " << (*trpart).numberOfHits() << " PSimHits"
0363               << " from " << (*trpart).g4Tracks().size() << " SimTrack:";
0364           for (TrackingParticle::g4t_iterator g4T = (*trpart).g4Track_begin(); g4T != (*trpart).g4Track_end(); ++g4T) {
0365             edm::LogVerbatim("MuonAssociatorByHitsHelper")
0366                 << "\t"
0367                 << " Id:" << (*g4T).trackId() << "/Evt:(" << (*g4T).eventId().event() << ","
0368                 << (*g4T).eventId().bunchCrossing() << ")";
0369           }
0370         } else {
0371           // print something only if this TrackingParticle shares some hits with
0372           // the current reco::Track
0373           if (global_nshared != 0)
0374             LogTrace("MuonAssociatorByHitsHelper")
0375                 << "\n\t"
0376                 << " NOT matched to TrackingParticle, with quality = " << global_quality
0377                 << " (tracker: " << tracker_quality << " / muon: " << muon_quality << ")"
0378                 << "\n"
0379                 << "   N shared hits = " << global_nshared << " (tracker: " << tracker_nshared
0380                 << " / muon: " << muon_nshared << ")";
0381         }
0382 
0383       }  //  loop over TrackingParticle
0384 
0385       if (!this_track_matched) {
0386         edm::LogVerbatim("MuonAssociatorByHitsHelper") << "\n"
0387                                                        << " NOT matched to any TrackingParticle";
0388       }
0389 
0390       edm::LogVerbatim("MuonAssociatorByHitsHelper")
0391           << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
0392              "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
0393           << "\n";
0394 
0395     }  //  if(n_matching_simhits != 0)
0396 
0397   }  // loop over reco::Track
0398 
0399   if (tC.empty())
0400     edm::LogVerbatim("MuonAssociatorByHitsHelper") << "0 reconstructed tracks (-->> 0 associated !)";
0401 
0402   for (IndexAssociation::iterator it = outputCollection.begin(), ed = outputCollection.end(); it != ed; ++it) {
0403     std::sort(it->second.begin(), it->second.end());
0404   }
0405   return outputCollection;
0406 }
0407 
0408 MuonAssociatorByHitsHelper::IndexAssociation MuonAssociatorByHitsHelper::associateSimToRecoIndices(
0409     const TrackHitsCollection &tC,
0410     const edm::RefVector<TrackingParticleCollection> &TPCollectionH,
0411     const Resources &resources) const {
0412   auto tTopo = resources.tTopo_;
0413   auto trackertruth = resources.trackerHitAssoc_;
0414   auto const &csctruth = *resources.cscHitAssoc_;
0415   auto const &dttruth = *resources.dtHitAssoc_;
0416   auto const &rpctruth = *resources.rpcHitAssoc_;
0417   auto const &gemtruth = *resources.gemHitAssoc_;
0418 
0419   int tracker_nshared = 0;
0420   int muon_nshared = 0;
0421   int global_nshared = 0;
0422 
0423   double tracker_quality = 0;
0424   double tracker_quality_cut;
0425   if (AbsoluteNumberOfHits_track)
0426     tracker_quality_cut = static_cast<double>(NHitCut_track);
0427   else
0428     tracker_quality_cut = EfficiencyCut_track;
0429 
0430   double muon_quality = 0;
0431   double muon_quality_cut;
0432   if (AbsoluteNumberOfHits_muon)
0433     muon_quality_cut = static_cast<double>(NHitCut_muon);
0434   else
0435     muon_quality_cut = EfficiencyCut_muon;
0436 
0437   double global_quality = 0;
0438 
0439   double tracker_purity = 0;
0440   double muon_purity = 0;
0441   double global_purity = 0;
0442 
0443   MapOfMatchedIds tracker_matchedIds_valid, muon_matchedIds_valid;
0444   MapOfMatchedIds tracker_matchedIds_INVALID, muon_matchedIds_INVALID;
0445 
0446   IndexAssociation outputCollection;
0447 
0448   bool printRtS = false;
0449 
0450   TrackingParticleCollection tPC;
0451   tPC.reserve(TPCollectionH.size());
0452   for (auto const &ref : TPCollectionH) {
0453     tPC.push_back(*ref);
0454   }
0455 
0456   bool any_trackingParticle_matched = false;
0457 
0458   int tindex = 0;
0459   for (TrackHitsCollection::const_iterator track = tC.begin(); track != tC.end(); track++, tindex++) {
0460     if (printRtS)
0461       edm::LogVerbatim("MuonAssociatorByHitsHelper")
0462           << "\n"
0463           << "reco::Track " << tindex << ", number of RecHits = " << (track->second - track->first) << "\n";
0464 
0465     tracker_matchedIds_valid.clear();
0466     muon_matchedIds_valid.clear();
0467 
0468     tracker_matchedIds_INVALID.clear();
0469     muon_matchedIds_INVALID.clear();
0470 
0471     int n_matching_simhits = 0;
0472 
0473     // all hits = valid +INVALID
0474     int n_all = 0;
0475     int n_tracker_all = 0;
0476     int n_dt_all = 0;
0477     int n_csc_all = 0;
0478     int n_rpc_all = 0;
0479     int n_gem_all = 0;
0480 
0481     int n_valid = 0;
0482     int n_tracker_valid = 0;
0483     int n_muon_valid = 0;
0484     int n_dt_valid = 0;
0485     int n_csc_valid = 0;
0486     int n_rpc_valid = 0;
0487     int n_gem_valid = 0;
0488 
0489     int n_tracker_matched_valid = 0;
0490     int n_muon_matched_valid = 0;
0491     int n_dt_matched_valid = 0;
0492     int n_csc_matched_valid = 0;
0493     int n_rpc_matched_valid = 0;
0494     int n_gem_matched_valid = 0;
0495 
0496     int n_INVALID = 0;
0497     int n_tracker_INVALID = 0;
0498     int n_muon_INVALID = 0;
0499     int n_dt_INVALID = 0;
0500     int n_csc_INVALID = 0;
0501     int n_rpc_INVALID = 0;
0502     int n_gem_INVALID = 0;
0503 
0504     int n_tracker_matched_INVALID = 0;
0505     int n_muon_matched_INVALID = 0;
0506     int n_dt_matched_INVALID = 0;
0507     int n_csc_matched_INVALID = 0;
0508     int n_rpc_matched_INVALID = 0;
0509     int n_gem_matched_INVALID = 0;
0510 
0511     getMatchedIds(tracker_matchedIds_valid,
0512                   muon_matchedIds_valid,
0513                   tracker_matchedIds_INVALID,
0514                   muon_matchedIds_INVALID,
0515                   n_tracker_valid,
0516                   n_dt_valid,
0517                   n_csc_valid,
0518                   n_rpc_valid,
0519                   n_gem_valid,
0520                   n_tracker_matched_valid,
0521                   n_dt_matched_valid,
0522                   n_csc_matched_valid,
0523                   n_rpc_matched_valid,
0524                   n_gem_matched_valid,
0525                   n_tracker_INVALID,
0526                   n_dt_INVALID,
0527                   n_csc_INVALID,
0528                   n_rpc_INVALID,
0529                   n_gem_INVALID,
0530                   n_tracker_matched_INVALID,
0531                   n_dt_matched_INVALID,
0532                   n_csc_matched_INVALID,
0533                   n_rpc_matched_INVALID,
0534                   n_gem_matched_INVALID,
0535                   track->first,
0536                   track->second,
0537                   trackertruth,
0538                   dttruth,
0539                   csctruth,
0540                   rpctruth,
0541                   gemtruth,
0542                   printRtS,
0543                   tTopo);
0544 
0545     n_matching_simhits = tracker_matchedIds_valid.size() + muon_matchedIds_valid.size() +
0546                          tracker_matchedIds_INVALID.size() + muon_matchedIds_INVALID.size();
0547 
0548     n_muon_valid = n_dt_valid + n_csc_valid + n_rpc_valid + n_gem_valid;
0549     n_valid = n_tracker_valid + n_muon_valid;
0550     n_muon_INVALID = n_dt_INVALID + n_csc_INVALID + n_rpc_INVALID + n_gem_INVALID;
0551     n_INVALID = n_tracker_INVALID + n_muon_INVALID;
0552 
0553     // all used hits (valid+INVALID), defined by UseTracker, UseMuon
0554     n_tracker_all = n_tracker_valid + n_tracker_INVALID;
0555     n_dt_all = n_dt_valid + n_dt_INVALID;
0556     n_csc_all = n_csc_valid + n_csc_INVALID;
0557     n_rpc_all = n_rpc_valid + n_rpc_INVALID;
0558     n_gem_all = n_gem_valid + n_gem_INVALID;
0559     n_all = n_valid + n_INVALID;
0560 
0561     n_muon_matched_valid = n_dt_matched_valid + n_csc_matched_valid + n_rpc_matched_valid + n_gem_matched_valid;
0562     n_muon_matched_INVALID =
0563         n_dt_matched_INVALID + n_csc_matched_INVALID + n_rpc_matched_INVALID + n_gem_matched_INVALID;
0564 
0565     // selected hits are set initially to valid hits
0566     int n_tracker_selected_hits = n_tracker_valid;
0567     int n_muon_selected_hits = n_muon_valid;
0568     int n_dt_selected_hits = n_dt_valid;
0569     int n_csc_selected_hits = n_csc_valid;
0570     int n_rpc_selected_hits = n_rpc_valid;
0571     int n_gem_selected_hits = n_gem_valid;
0572 
0573     // matched hits are a subsample of the selected hits
0574     int n_tracker_matched = n_tracker_matched_valid;
0575     int n_muon_matched = n_muon_matched_valid;
0576     int n_dt_matched = n_dt_matched_valid;
0577     int n_csc_matched = n_csc_matched_valid;
0578     int n_rpc_matched = n_rpc_matched_valid;
0579     int n_gem_matched = n_gem_matched_valid;
0580 
0581     std::string InvMuonHits, ZeroHitMuon;
0582 
0583     if (includeZeroHitMuons && n_muon_valid == 0 && n_muon_INVALID != 0) {
0584       // selected muon hits = INVALID when (useZeroHitMuons == True) and track
0585       // has no valid muon hits
0586 
0587       InvMuonHits = " ***INVALID MUON HITS***";
0588       ZeroHitMuon = " ***ZERO-HIT MUON***";
0589 
0590       n_muon_selected_hits = n_muon_INVALID;
0591       n_dt_selected_hits = n_dt_INVALID;
0592       n_csc_selected_hits = n_csc_INVALID;
0593       n_rpc_selected_hits = n_rpc_INVALID;
0594       n_gem_selected_hits = n_gem_INVALID;
0595 
0596       n_muon_matched = n_muon_matched_INVALID;
0597       n_dt_matched = n_dt_matched_INVALID;
0598       n_csc_matched = n_csc_matched_INVALID;
0599       n_rpc_matched = n_rpc_matched_INVALID;
0600       n_gem_matched = n_gem_matched_INVALID;
0601     }
0602 
0603     int n_selected_hits = n_tracker_selected_hits + n_muon_selected_hits;
0604     int n_matched = n_tracker_matched + n_muon_matched;
0605 
0606     if (printRtS)
0607       edm::LogVerbatim("MuonAssociatorByHitsHelper")
0608           << "\n"
0609           << "# TrackingRecHits: " << (track->second - track->first) << "\n"
0610           << "# used RecHits     = " << n_all << " (" << n_tracker_all << "/" << n_dt_all << "/" << n_csc_all << "/"
0611           << n_rpc_all << "/" << n_gem_all << " in Tracker/DT/CSC/RPC/GEM)"
0612           << ", obtained from " << n_matching_simhits << " SimHits"
0613           << "\n"
0614           << "# selected RecHits = " << n_selected_hits << " (" << n_tracker_selected_hits << "/" << n_dt_selected_hits
0615           << "/" << n_csc_selected_hits << "/" << n_rpc_selected_hits << "/" << n_gem_selected_hits
0616           << " in Tracker/DT/CSC/RPC/GEM)" << InvMuonHits << "\n"
0617           << "# matched RecHits = " << n_matched << " (" << n_tracker_matched << "/" << n_dt_matched << "/"
0618           << n_csc_matched << "/" << n_rpc_matched << "/" << n_gem_matched << " in Tracker/DT/CSC/RPC/GEM)";
0619 
0620     if (printRtS && n_all > 0 && n_matching_simhits == 0)
0621       edm::LogVerbatim("MuonAssociatorByHitsHelper")
0622           << "*** WARNING in MuonAssociatorByHitsHelper::associateSimToReco: "
0623              "no matching PSimHit found for this reco::Track !";
0624 
0625     if (n_matching_simhits != 0) {
0626       int tpindex = 0;
0627       for (TrackingParticleCollection::const_iterator trpart = tPC.begin(); trpart != tPC.end(); ++trpart, ++tpindex) {
0628         //  int n_tracker_simhits = 0;
0629         int n_tracker_recounted_simhits = 0;
0630         int n_muon_simhits = 0;
0631         int n_global_simhits = 0;
0632 
0633         int n_tracker_selected_simhits = 0;
0634         int n_muon_selected_simhits = 0;
0635         int n_global_selected_simhits = 0;
0636 
0637         // shared hits are counted over the selected ones
0638         tracker_nshared = getShared(tracker_matchedIds_valid, trpart);
0639         muon_nshared = getShared(muon_matchedIds_valid, trpart);
0640 
0641         if (includeZeroHitMuons && n_muon_valid == 0 && n_muon_INVALID != 0)
0642           muon_nshared = getShared(muon_matchedIds_INVALID, trpart);
0643 
0644         global_nshared = tracker_nshared + muon_nshared;
0645         if (global_nshared == 0)
0646           continue;  // if this TP shares no hits with the current reco::Track
0647                      // loop over
0648 
0649         // adapt to new TP interface: this gives the total number of hits in
0650         // tracker
0651         //   should reproduce the behaviour of UseGrouped=UseSplitting=.true.
0652         n_tracker_recounted_simhits = trpart->numberOfTrackerHits();
0653         //   numberOfHits() gives the total number of hits (tracker + muons)
0654         n_muon_simhits = trpart->numberOfHits() - trpart->numberOfTrackerHits();
0655 
0656         // Handle the case of TrackingParticles that don't have PSimHits inside,
0657         // e.g. because they were made on RECOSIM only.
0658         if (trpart->numberOfHits() == 0) {
0659           // FIXME this can be made better, counting the digiSimLinks associated
0660           // to this TP, but perhaps it's not worth it
0661           // n_tracker_recounted_simhits = tracker_nshared;
0662           // n_muon_simhits = muon_nshared;
0663           // ---> on RECOSIM when the AbsoluteNumberOfHits_muon=True this always
0664           // obtains quality=1,
0665           //      hence no sorting is possible in case of duplicate matchings
0666           // ---> reset these variables to 1 so to keep the number of shared
0667           // hits as ranking criterion
0668           n_tracker_recounted_simhits = 1;
0669           n_muon_simhits = 1;
0670         }
0671         n_global_simhits = n_tracker_recounted_simhits + n_muon_simhits;
0672 
0673         if (UseMuon) {
0674           n_muon_selected_simhits = n_muon_simhits;
0675           n_global_selected_simhits = n_muon_selected_simhits;
0676         }
0677         if (UseTracker) {
0678           n_tracker_selected_simhits = n_tracker_recounted_simhits;
0679           n_global_selected_simhits += n_tracker_selected_simhits;
0680         }
0681 
0682         if (AbsoluteNumberOfHits_track)
0683           tracker_quality = static_cast<double>(tracker_nshared);
0684         else if (n_tracker_selected_simhits != 0)
0685           tracker_quality = static_cast<double>(tracker_nshared) / static_cast<double>(n_tracker_selected_simhits);
0686         else
0687           tracker_quality = 0;
0688 
0689         if (AbsoluteNumberOfHits_muon)
0690           muon_quality = static_cast<double>(muon_nshared);
0691         else if (n_muon_selected_simhits != 0)
0692           muon_quality = static_cast<double>(muon_nshared) / static_cast<double>(n_muon_selected_simhits);
0693         else
0694           muon_quality = 0;
0695 
0696         // global_quality used to order the matching tracks
0697         if (n_global_selected_simhits != 0) {
0698           if (AbsoluteNumberOfHits_muon && AbsoluteNumberOfHits_track)
0699             global_quality = global_nshared;
0700           else
0701             global_quality = static_cast<double>(global_nshared) / static_cast<double>(n_global_selected_simhits);
0702         } else
0703           global_quality = 0;
0704 
0705         // global purity
0706         if (n_selected_hits != 0) {
0707           if (AbsoluteNumberOfHits_muon && AbsoluteNumberOfHits_track)
0708             global_purity = global_nshared;
0709           else
0710             global_purity = static_cast<double>(global_nshared) / static_cast<double>(n_selected_hits);
0711         } else
0712           global_purity = 0;
0713 
0714         bool trackerOk = false;
0715         if (n_tracker_selected_hits != 0) {
0716           if (tracker_quality > tracker_quality_cut)
0717             trackerOk = true;
0718 
0719           tracker_purity = static_cast<double>(tracker_nshared) / static_cast<double>(n_tracker_selected_hits);
0720           if (AbsoluteNumberOfHits_track)
0721             tracker_purity = static_cast<double>(tracker_nshared);
0722 
0723           if ((!AbsoluteNumberOfHits_track) && tracker_purity <= PurityCut_track)
0724             trackerOk = false;
0725 
0726           // if a track has just 3 hits in the Tracker we require that all 3
0727           // hits are shared
0728           if (ThreeHitTracksAreSpecial && n_tracker_selected_hits == 3 && tracker_nshared < 3)
0729             trackerOk = false;
0730         }
0731 
0732         bool muonOk = false;
0733         if (n_muon_selected_hits != 0) {
0734           if (muon_quality > muon_quality_cut)
0735             muonOk = true;
0736 
0737           muon_purity = static_cast<double>(muon_nshared) / static_cast<double>(n_muon_selected_hits);
0738           if (AbsoluteNumberOfHits_muon)
0739             muon_purity = static_cast<double>(muon_nshared);
0740 
0741           if ((!AbsoluteNumberOfHits_muon) && muon_purity <= PurityCut_muon)
0742             muonOk = false;
0743         }
0744 
0745         // (matchOk) has to account for different track types (tracker-only,
0746         // standalone muons, global muons)
0747         bool matchOk = trackerOk || muonOk;
0748 
0749         // only for global tracks: match both tracker and muon stub (if acceptOneStubMatchings==false)
0750         // depending on the muon selection reject tracks with only one stub (if rejectBadGlobal==true)
0751         //
0752         if (UseTracker && UseMuon && !acceptOneStubMatchings &&
0753             ((n_tracker_selected_hits != 0 && n_muon_selected_hits != 0) || rejectBadGlobal))
0754           matchOk = trackerOk && muonOk;
0755 
0756         if (matchOk) {
0757           outputCollection[tpindex].push_back(IndexMatch(tindex, global_quality));
0758           any_trackingParticle_matched = true;
0759 
0760           edm::LogVerbatim("MuonAssociatorByHitsHelper")
0761               << "*************************************************************"
0762                  "***********************************************************"
0763               << "\n"
0764               << "TrackingParticle:  q = " << (*trpart).charge() << ", p = " << (*trpart).p()
0765               << ", pT = " << (*trpart).pt() << ", eta = " << (*trpart).eta() << ", phi = " << (*trpart).phi() << "\n"
0766               << " pdg code = " << (*trpart).pdgId() << ", made of " << (*trpart).numberOfHits()
0767               << " PSimHits, recounted " << n_global_simhits << " PSimHits"
0768               << " (tracker:" << n_tracker_recounted_simhits << "/muons:" << n_muon_simhits << ")"
0769               << ", from " << (*trpart).g4Tracks().size() << " SimTrack:";
0770           for (TrackingParticle::g4t_iterator g4T = (*trpart).g4Track_begin(); g4T != (*trpart).g4Track_end(); ++g4T) {
0771             edm::LogVerbatim("MuonAssociatorByHitsHelper")
0772                 << " Id:" << (*g4T).trackId() << "/Evt:(" << (*g4T).eventId().event() << ","
0773                 << (*g4T).eventId().bunchCrossing() << ")";
0774           }
0775           edm::LogVerbatim("MuonAssociatorByHitsHelper")
0776               << "\t selected " << n_global_selected_simhits << " PSimHits"
0777               << " (tracker:" << n_tracker_selected_simhits << "/muons:" << n_muon_selected_simhits << ")"
0778               << "\n\t **MATCHED** with quality = " << global_quality << " (tracker: " << tracker_quality
0779               << " / muon: " << muon_quality << ")"
0780               << "\n\t               and purity = " << global_purity << " (tracker: " << tracker_purity
0781               << " / muon: " << muon_purity << ")"
0782               << "\n\t   N shared hits = " << global_nshared << " (tracker: " << tracker_nshared
0783               << " / muon: " << muon_nshared << ")"
0784               << "\n"
0785               << "   to: reco::Track " << tindex << ZeroHitMuon << "\n\t"
0786               << " made of " << n_selected_hits << " RecHits (tracker:" << n_tracker_valid
0787               << "/muons:" << n_muon_selected_hits << ")";
0788         } else {
0789           // print something only if this TrackingParticle shares some hits with
0790           // the current reco::Track
0791           if (global_nshared != 0) {
0792             if (printRtS)
0793               edm::LogVerbatim("MuonAssociatorByHitsHelper")
0794                   << "*********************************************************"
0795                      "*********************************************************"
0796                      "******"
0797                   << "\n"
0798                   << "TrackingParticle:  q = " << (*trpart).charge() << ", p = " << (*trpart).p()
0799                   << ", pT = " << (*trpart).pt() << ", eta = " << (*trpart).eta() << ", phi = " << (*trpart).phi()
0800                   << "\n"
0801                   << " pdg code = " << (*trpart).pdgId() << ", made of " << (*trpart).numberOfHits()
0802                   << " PSimHits, recounted " << n_global_simhits << " PSimHits"
0803                   << " (tracker:" << n_tracker_recounted_simhits << "/muons:" << n_muon_simhits << ")"
0804                   << ", from " << (*trpart).g4Tracks().size() << " SimTrack:";
0805             for (TrackingParticle::g4t_iterator g4T = (*trpart).g4Track_begin(); g4T != (*trpart).g4Track_end();
0806                  ++g4T) {
0807               if (printRtS)
0808                 edm::LogVerbatim("MuonAssociatorByHitsHelper")
0809                     << " Id:" << (*g4T).trackId() << "/Evt:(" << (*g4T).eventId().event() << ","
0810                     << (*g4T).eventId().bunchCrossing() << ")";
0811             }
0812             if (printRtS)
0813               edm::LogVerbatim("MuonAssociatorByHitsHelper")
0814                   << "\t selected " << n_global_selected_simhits << " PSimHits"
0815                   << " (tracker:" << n_tracker_selected_simhits << "/muons:" << n_muon_selected_simhits << ")"
0816                   << "\n\t NOT matched  to reco::Track " << tindex << ZeroHitMuon
0817                   << " with quality = " << global_quality << " (tracker: " << tracker_quality
0818                   << " / muon: " << muon_quality << ")"
0819                   << "\n\t and purity = " << global_purity << " (tracker: " << tracker_purity
0820                   << " / muon: " << muon_purity << ")"
0821                   << "\n\t     N shared hits = " << global_nshared << " (tracker: " << tracker_nshared
0822                   << " / muon: " << muon_nshared << ")";
0823           }
0824         }
0825       }  // loop over TrackingParticle's
0826     }  // if(n_matching_simhits != 0)
0827   }  // loop over reco Tracks
0828 
0829   if (!any_trackingParticle_matched) {
0830     edm::LogVerbatim("MuonAssociatorByHitsHelper")
0831         << "\n"
0832         << "*******************************************************************"
0833            "*****************************************************"
0834         << "\n NO TrackingParticle associated to ANY input reco::Track ! \n"
0835         << "*******************************************************************"
0836            "*****************************************************"
0837         << "\n";
0838   } else {
0839     edm::LogVerbatim("MuonAssociatorByHitsHelper")
0840         << "*******************************************************************"
0841            "*****************************************************"
0842         << "\n";
0843   }
0844 
0845   for (IndexAssociation::iterator it = outputCollection.begin(), ed = outputCollection.end(); it != ed; ++it) {
0846     std::sort(it->second.begin(), it->second.end());
0847   }
0848   return outputCollection;
0849 }
0850 
0851 void MuonAssociatorByHitsHelper::getMatchedIds(MapOfMatchedIds &tracker_matchedIds_valid,
0852                                                MapOfMatchedIds &muon_matchedIds_valid,
0853                                                MapOfMatchedIds &tracker_matchedIds_INVALID,
0854                                                MapOfMatchedIds &muon_matchedIds_INVALID,
0855                                                int &n_tracker_valid,
0856                                                int &n_dt_valid,
0857                                                int &n_csc_valid,
0858                                                int &n_rpc_valid,
0859                                                int &n_gem_valid,
0860                                                int &n_tracker_matched_valid,
0861                                                int &n_dt_matched_valid,
0862                                                int &n_csc_matched_valid,
0863                                                int &n_rpc_matched_valid,
0864                                                int &n_gem_matched_valid,
0865                                                int &n_tracker_INVALID,
0866                                                int &n_dt_INVALID,
0867                                                int &n_csc_INVALID,
0868                                                int &n_rpc_INVALID,
0869                                                int &n_gem_INVALID,
0870                                                int &n_tracker_matched_INVALID,
0871                                                int &n_dt_matched_INVALID,
0872                                                int &n_csc_matched_INVALID,
0873                                                int &n_rpc_matched_INVALID,
0874                                                int &n_gem_matched_INVALID,
0875                                                trackingRecHit_iterator begin,
0876                                                trackingRecHit_iterator end,
0877                                                const TrackerHitAssociator *trackertruth,
0878                                                const DTHitAssociator &dttruth,
0879                                                const CSCHitAssociator &csctruth,
0880                                                const RPCHitAssociator &rpctruth,
0881                                                const GEMHitAssociator &gemtruth,
0882                                                bool printRtS,
0883                                                const TrackerTopology *tTopo) const
0884 
0885 {
0886   tracker_matchedIds_valid.clear();
0887   muon_matchedIds_valid.clear();
0888 
0889   tracker_matchedIds_INVALID.clear();
0890   muon_matchedIds_INVALID.clear();
0891 
0892   n_tracker_valid = 0;
0893   n_dt_valid = 0;
0894   n_csc_valid = 0;
0895   n_rpc_valid = 0;
0896   n_gem_valid = 0;
0897 
0898   n_tracker_matched_valid = 0;
0899   n_dt_matched_valid = 0;
0900   n_csc_matched_valid = 0;
0901   n_rpc_matched_valid = 0;
0902   n_gem_matched_valid = 0;
0903 
0904   n_tracker_INVALID = 0;
0905   n_dt_INVALID = 0;
0906   n_csc_INVALID = 0;
0907   n_rpc_INVALID = 0;
0908   n_gem_INVALID = 0;
0909 
0910   n_tracker_matched_INVALID = 0;
0911   n_dt_matched_INVALID = 0;
0912   n_csc_matched_INVALID = 0;
0913   n_rpc_matched_INVALID = 0;
0914   n_gem_matched_INVALID = 0;
0915 
0916   std::vector<SimHitIdpr> SimTrackIds;
0917 
0918   // main loop on TrackingRecHits
0919   int iloop = 0;
0920   int iH = -1;
0921   for (trackingRecHit_iterator it = begin; it != end; it++, iloop++) {
0922     stringstream hit_index;
0923     hit_index << iloop;
0924 
0925     const TrackingRecHit *hitp = getHitPtr(it);
0926     DetId geoid = hitp->geographicalId();
0927 
0928     unsigned int detid = geoid.rawId();
0929     stringstream detector_id;
0930     detector_id << detid;
0931 
0932     string hitlog = "TrackingRecHit " + hit_index.str();
0933     string wireidlog;
0934     std::vector<string> DTSimHits;
0935 
0936     DetId::Detector det = geoid.det();
0937     int subdet = geoid.subdetId();
0938 
0939     bool valid_Hit = hitp->isValid();
0940 
0941     // Si-Tracker Hits
0942     if (det == DetId::Tracker && UseTracker) {
0943       stringstream detector_id;
0944       detector_id << tTopo->print(detid);
0945 
0946       if (valid_Hit)
0947         hitlog = hitlog + " -Tracker - detID = " + detector_id.str();
0948       else
0949         hitlog = hitlog + " *** INVALID ***" + " -Tracker - detID = " + detector_id.str();
0950 
0951       iH++;
0952       SimTrackIds = trackertruth->associateHitId(*hitp);
0953 
0954       if (valid_Hit) {
0955         n_tracker_valid++;
0956 
0957         if (!SimTrackIds.empty()) {
0958           n_tracker_matched_valid++;
0959           // tracker_matchedIds_valid[iH] = SimTrackIds;
0960           tracker_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
0961         }
0962       } else {
0963         n_tracker_INVALID++;
0964 
0965         if (!SimTrackIds.empty()) {
0966           n_tracker_matched_INVALID++;
0967           // tracker_matchedIds_INVALID[iH] = SimTrackIds;
0968           tracker_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
0969         }
0970       }
0971     }
0972     // Muon detector Hits
0973     else if (det == DetId::Muon && UseMuon) {
0974       // DT Hits
0975       if (subdet == MuonSubdetId::DT) {
0976         DTWireId dtdetid = DTWireId(detid);
0977         stringstream dt_detector_id;
0978         dt_detector_id << dtdetid;
0979         if (valid_Hit)
0980           hitlog = hitlog + " -Muon DT - detID = " + dt_detector_id.str();
0981         else
0982           hitlog = hitlog + " *** INVALID ***" + " -Muon DT - detID = " + dt_detector_id.str();
0983 
0984         const DTRecHit1D *dtrechit = dynamic_cast<const DTRecHit1D *>(hitp);
0985 
0986         // single DT hits
0987         if (dtrechit) {
0988           iH++;
0989           SimTrackIds = dttruth.associateDTHitId(dtrechit);
0990 
0991           if (valid_Hit) {
0992             n_dt_valid++;
0993 
0994             if (!SimTrackIds.empty()) {
0995               n_dt_matched_valid++;
0996               // muon_matchedIds_valid[iH] = SimTrackIds;
0997               muon_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
0998             }
0999           } else {
1000             n_dt_INVALID++;
1001 
1002             if (!SimTrackIds.empty()) {
1003               n_dt_matched_INVALID++;
1004               // muon_matchedIds_INVALID[iH] = SimTrackIds;
1005               muon_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
1006             }
1007           }
1008 
1009           if (dumpDT) {
1010             DTWireId wireid = dtrechit->wireId();
1011             stringstream wid;
1012             wid << wireid;
1013             std::vector<PSimHit> dtSimHits = dttruth.associateHit(*hitp);
1014 
1015             stringstream ndthits;
1016             ndthits << dtSimHits.size();
1017             wireidlog = "\t DTWireId :" + wid.str() + ", " + ndthits.str() + " associated PSimHit :";
1018 
1019             for (unsigned int j = 0; j < dtSimHits.size(); j++) {
1020               stringstream index;
1021               index << j;
1022               stringstream simhit;
1023               simhit << dtSimHits[j];
1024               string simhitlog = "\t\t PSimHit " + index.str() + ": " + simhit.str();
1025               DTSimHits.push_back(simhitlog);
1026             }
1027           }  // if (dumpDT)
1028         }
1029 
1030         // DT segments
1031         else {
1032           const DTRecSegment4D *dtsegment = dynamic_cast<const DTRecSegment4D *>(hitp);
1033 
1034           if (dtsegment) {
1035             std::vector<const TrackingRecHit *> componentHits, phiHits, zHits;
1036             if (dtsegment->hasPhi()) {
1037               phiHits = dtsegment->phiSegment()->recHits();
1038               componentHits.insert(componentHits.end(), phiHits.begin(), phiHits.end());
1039             }
1040             if (dtsegment->hasZed()) {
1041               zHits = dtsegment->zSegment()->recHits();
1042               componentHits.insert(componentHits.end(), zHits.begin(), zHits.end());
1043             }
1044             if (printRtS)
1045               edm::LogVerbatim("MuonAssociatorByHitsHelper")
1046                   << "\n\t this TrackingRecHit is a DTRecSegment4D with " << componentHits.size()
1047                   << " hits (phi:" << phiHits.size() << ", z:" << zHits.size() << ")";
1048 
1049             SimTrackIds.clear();
1050             std::vector<SimHitIdpr> i_SimTrackIds;
1051             int i_compHit = 0;
1052             for (std::vector<const TrackingRecHit *>::const_iterator ithit = componentHits.begin();
1053                  ithit != componentHits.end();
1054                  ++ithit) {
1055               i_compHit++;
1056 
1057               const DTRecHit1D *dtrechit1D = dynamic_cast<const DTRecHit1D *>(*ithit);
1058 
1059               i_SimTrackIds.clear();
1060               if (dtrechit1D) {
1061                 iH++;
1062                 i_SimTrackIds = dttruth.associateDTHitId(dtrechit1D);
1063 
1064                 if (valid_Hit) {
1065                   // validity check is on the segment, but hits are counted
1066                   // one-by-one
1067                   n_dt_valid++;
1068 
1069                   if (!i_SimTrackIds.empty()) {
1070                     n_dt_matched_valid++;
1071                     // muon_matchedIds_valid[iH] = i_SimTrackIds;
1072                     muon_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, i_SimTrackIds));
1073                   }
1074                 } else {
1075                   n_dt_INVALID++;
1076 
1077                   if (!i_SimTrackIds.empty()) {
1078                     n_dt_matched_INVALID++;
1079                     // muon_matchedIds_INVALID[iH] = i_SimTrackIds;
1080                     muon_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, i_SimTrackIds));
1081                   }
1082                 }
1083               } else if (printRtS)
1084                 edm::LogVerbatim("MuonAssociatorByHitsHelper") << "*** WARNING in "
1085                                                                   "MuonAssociatorByHitsHelper::getMatchedIds, null "
1086                                                                   "dynamic_cast of a DT TrackingRecHit !";
1087 
1088               unsigned int i_detid = (*ithit)->geographicalId().rawId();
1089               DTWireId i_dtdetid = DTWireId(i_detid);
1090 
1091               stringstream i_dt_detector_id;
1092               i_dt_detector_id << i_dtdetid;
1093 
1094               stringstream i_ss;
1095               i_ss << "\t\t hit " << i_compHit << " -Muon DT - detID = " << i_dt_detector_id.str();
1096 
1097               string i_hitlog = i_ss.str();
1098               i_hitlog = i_hitlog + write_matched_simtracks(i_SimTrackIds);
1099               if (printRtS)
1100                 edm::LogVerbatim("MuonAssociatorByHitsHelper") << i_hitlog;
1101 
1102               SimTrackIds.insert(SimTrackIds.end(), i_SimTrackIds.begin(), i_SimTrackIds.end());
1103             }
1104           }  // if (dtsegment)
1105 
1106           else if (printRtS)
1107             edm::LogVerbatim("MuonAssociatorByHitsHelper")
1108                 << "*** WARNING in MuonAssociatorByHitsHelper::getMatchedIds, "
1109                    "DT TrackingRecHit is neither DTRecHit1D nor DTRecSegment4D "
1110                    "! ";
1111         }
1112       }
1113 
1114       // CSC Hits
1115       else if (subdet == MuonSubdetId::CSC) {
1116         CSCDetId cscdetid = CSCDetId(detid);
1117         stringstream csc_detector_id;
1118         csc_detector_id << cscdetid;
1119         if (valid_Hit)
1120           hitlog = hitlog + " -Muon CSC- detID = " + csc_detector_id.str();
1121         else
1122           hitlog = hitlog + " *** INVALID ***" + " -Muon CSC- detID = " + csc_detector_id.str();
1123 
1124         const CSCRecHit2D *cscrechit = dynamic_cast<const CSCRecHit2D *>(hitp);
1125 
1126         // single CSC hits
1127         if (cscrechit) {
1128           iH++;
1129           SimTrackIds = csctruth.associateCSCHitId(cscrechit);
1130 
1131           if (valid_Hit) {
1132             n_csc_valid++;
1133 
1134             if (!SimTrackIds.empty()) {
1135               n_csc_matched_valid++;
1136               // muon_matchedIds_valid[iH] = SimTrackIds;
1137               muon_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
1138             }
1139           } else {
1140             n_csc_INVALID++;
1141 
1142             if (!SimTrackIds.empty()) {
1143               n_csc_matched_INVALID++;
1144               // muon_matchedIds_INVALID[iH] = SimTrackIds;
1145               muon_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
1146             }
1147           }
1148         }
1149 
1150         // CSC segments
1151         else {
1152           const CSCSegment *cscsegment = dynamic_cast<const CSCSegment *>(hitp);
1153 
1154           if (cscsegment) {
1155             std::vector<const TrackingRecHit *> componentHits = cscsegment->recHits();
1156             if (printRtS)
1157               edm::LogVerbatim("MuonAssociatorByHitsHelper")
1158                   << "\n\t this TrackingRecHit is a CSCSegment with " << componentHits.size() << " hits";
1159 
1160             SimTrackIds.clear();
1161             std::vector<SimHitIdpr> i_SimTrackIds;
1162             int i_compHit = 0;
1163             for (std::vector<const TrackingRecHit *>::const_iterator ithit = componentHits.begin();
1164                  ithit != componentHits.end();
1165                  ++ithit) {
1166               i_compHit++;
1167 
1168               const CSCRecHit2D *cscrechit2D = dynamic_cast<const CSCRecHit2D *>(*ithit);
1169 
1170               i_SimTrackIds.clear();
1171               if (cscrechit2D) {
1172                 iH++;
1173                 i_SimTrackIds = csctruth.associateCSCHitId(cscrechit2D);
1174 
1175                 if (valid_Hit) {
1176                   // validity check is on the segment, but hits are counted
1177                   // one-by-one
1178                   n_csc_valid++;
1179 
1180                   if (!i_SimTrackIds.empty()) {
1181                     n_csc_matched_valid++;
1182                     // muon_matchedIds_valid[iH] =  i_SimTrackIds;
1183                     muon_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, i_SimTrackIds));
1184                   }
1185                 } else {
1186                   n_csc_INVALID++;
1187 
1188                   if (!i_SimTrackIds.empty()) {
1189                     n_csc_matched_INVALID++;
1190                     // muon_matchedIds_INVALID[iH] =  i_SimTrackIds;
1191                     muon_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, i_SimTrackIds));
1192                   }
1193                 }
1194               } else if (printRtS)
1195                 edm::LogVerbatim("MuonAssociatorByHitsHelper") << "*** WARNING in "
1196                                                                   "MuonAssociatorByHitsHelper::getMatchedIds, null "
1197                                                                   "dynamic_cast of a CSC TrackingRecHit !";
1198 
1199               unsigned int i_detid = (*ithit)->geographicalId().rawId();
1200               CSCDetId i_cscdetid = CSCDetId(i_detid);
1201 
1202               stringstream i_csc_detector_id;
1203               i_csc_detector_id << i_cscdetid;
1204 
1205               stringstream i_ss;
1206               i_ss << "\t\t hit " << i_compHit << " -Muon CSC- detID = " << i_csc_detector_id.str();
1207 
1208               string i_hitlog = i_ss.str();
1209               i_hitlog = i_hitlog + write_matched_simtracks(i_SimTrackIds);
1210               if (printRtS)
1211                 edm::LogVerbatim("MuonAssociatorByHitsHelper") << i_hitlog;
1212 
1213               SimTrackIds.insert(SimTrackIds.end(), i_SimTrackIds.begin(), i_SimTrackIds.end());
1214             }
1215           }  // if (cscsegment)
1216 
1217           else if (printRtS)
1218             edm::LogVerbatim("MuonAssociatorByHitsHelper")
1219                 << "*** WARNING in MuonAssociatorByHitsHelper::getMatchedIds, "
1220                    "CSC TrackingRecHit is neither CSCRecHit2D nor CSCSegment "
1221                    "! ";
1222         }
1223       }
1224 
1225       // RPC Hits
1226       else if (subdet == MuonSubdetId::RPC) {
1227         RPCDetId rpcdetid = RPCDetId(detid);
1228         stringstream rpc_detector_id;
1229         rpc_detector_id << rpcdetid;
1230         if (valid_Hit)
1231           hitlog = hitlog + " -Muon RPC- detID = " + rpc_detector_id.str();
1232         else
1233           hitlog = hitlog + " *** INVALID ***" + " -Muon RPC- detID = " + rpc_detector_id.str();
1234 
1235         iH++;
1236         SimTrackIds = rpctruth.associateRecHit(*hitp);
1237 
1238         if (valid_Hit) {
1239           n_rpc_valid++;
1240 
1241           if (!SimTrackIds.empty()) {
1242             n_rpc_matched_valid++;
1243             // muon_matchedIds_valid[iH] = SimTrackIds;
1244             muon_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
1245           }
1246         } else {
1247           n_rpc_INVALID++;
1248 
1249           if (!SimTrackIds.empty()) {
1250             n_rpc_matched_INVALID++;
1251             // muon_matchedIds_INVALID[iH] = SimTrackIds;
1252             muon_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
1253           }
1254         }
1255       }
1256 
1257       // GEM Hits
1258       else if (subdet == MuonSubdetId::GEM) {
1259         GEMDetId gemdetid = GEMDetId(detid);
1260         stringstream gem_detector_id;
1261         gem_detector_id << gemdetid;
1262         if (valid_Hit)
1263           hitlog = hitlog + " -Muon GEM- detID = " + gem_detector_id.str();
1264         else
1265           hitlog = hitlog + " *** INVALID ***" + " -Muon GEM- detID = " + gem_detector_id.str();
1266 
1267         const GEMRecHit *gemrechit = dynamic_cast<const GEMRecHit *>(hitp);
1268         if (gemrechit) {
1269           iH++;
1270           SimTrackIds = gemtruth.associateRecHit(gemrechit);
1271 
1272           if (valid_Hit) {
1273             n_gem_valid++;
1274 
1275             if (!SimTrackIds.empty()) {
1276               n_gem_matched_valid++;
1277               // muon_matchedIds_valid[iH] = SimTrackIds;
1278               muon_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
1279             }
1280           } else {
1281             n_gem_INVALID++;
1282 
1283             if (!SimTrackIds.empty()) {
1284               n_gem_matched_INVALID++;
1285               // muon_matchedIds_INVALID[iH] = SimTrackIds;
1286               muon_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, SimTrackIds));
1287             }
1288           }
1289         } else {
1290           const GEMSegment *gemsegment = dynamic_cast<const GEMSegment *>(hitp);
1291           if (gemsegment) {
1292             std::vector<const TrackingRecHit *> componentHits = gemsegment->recHits();
1293             if (printRtS)
1294               edm::LogVerbatim("MuonAssociatorByHitsHelper")
1295                   << "\n\t this TrackingRecHit is a GEMSegment with " << componentHits.size() << " hits";
1296 
1297             SimTrackIds.clear();
1298             std::vector<SimHitIdpr> i_SimTrackIds;
1299 
1300             for (auto const &ithit : componentHits) {
1301               const GEMRecHit *gemrechitseg = dynamic_cast<const GEMRecHit *>(ithit);
1302 
1303               i_SimTrackIds.clear();
1304               if (gemrechitseg) {
1305                 iH++;
1306                 i_SimTrackIds = gemtruth.associateRecHit(gemrechitseg);
1307 
1308                 if (valid_Hit) {
1309                   // validity check is on the segment, but hits are counted
1310                   // one-by-one
1311                   n_gem_valid++;
1312 
1313                   if (!i_SimTrackIds.empty()) {
1314                     n_gem_matched_valid++;
1315                     // muon_matchedIds_valid[iH] =  i_SimTrackIds;
1316                     muon_matchedIds_valid.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, i_SimTrackIds));
1317                   }
1318                 } else {
1319                   n_gem_INVALID++;
1320 
1321                   if (!i_SimTrackIds.empty()) {
1322                     n_gem_matched_INVALID++;
1323                     // muon_matchedIds_INVALID[iH] =  i_SimTrackIds;
1324                     muon_matchedIds_INVALID.emplace_back(std::make_unique<uint_SimHitIdpr_pair>(iH, i_SimTrackIds));
1325                   }
1326                 }
1327               } else if (printRtS)
1328                 edm::LogVerbatim("MuonAssociatorByHitsHelper") << "*** WARNING in "
1329                                                                   "MuonAssociatorByHitsHelper::getMatchedIds, null "
1330                                                                   "dynamic_cast of a GEM TrackingRecHit !";
1331 
1332               if (printRtS) {
1333                 unsigned int i_detid = ithit->geographicalId().rawId();
1334                 GEMDetId i_gemdetid = GEMDetId(i_detid);
1335 
1336                 string i_hitlog = std::to_string(i_gemdetid);
1337                 i_hitlog = i_hitlog + write_matched_simtracks(i_SimTrackIds);
1338                 edm::LogVerbatim("MuonAssociatorByHitsHelper") << i_hitlog;
1339               }
1340 
1341               SimTrackIds.insert(SimTrackIds.end(), i_SimTrackIds.begin(), i_SimTrackIds.end());
1342             }
1343           }  // if (gemsegment)
1344 
1345         }  // if (gemrechit
1346       } else if (printRtS)
1347         edm::LogVerbatim("MuonAssociatorByHitsHelper")
1348             << "TrackingRecHit " << iloop << "  *** WARNING *** Unexpected Hit from Detector = " << det;
1349     }  // end if (det == DetId::Muon && UseMuon)
1350     else
1351       continue;
1352 
1353     hitlog = hitlog + write_matched_simtracks(SimTrackIds);
1354 
1355     if (printRtS)
1356       edm::LogVerbatim("MuonAssociatorByHitsHelper") << hitlog;
1357     if (printRtS && dumpDT && det == DetId::Muon && subdet == MuonSubdetId::DT) {
1358       edm::LogVerbatim("MuonAssociatorByHitsHelper") << wireidlog;
1359       for (unsigned int j = 0; j < DTSimHits.size(); j++) {
1360         edm::LogVerbatim("MuonAssociatorByHitsHelper") << DTSimHits[j];
1361       }
1362     }
1363 
1364   }  // trackingRecHit loop
1365 }
1366 
1367 int MuonAssociatorByHitsHelper::getShared(MapOfMatchedIds &matchedIds,
1368                                           TrackingParticleCollection::const_iterator trpart) const {
1369   int nshared = 0;
1370   const std::vector<SimTrack> &g4Tracks = trpart->g4Tracks();
1371 
1372   // map is indexed over the rechits of the reco::Track (no double-countings
1373   // allowed)
1374   for (MapOfMatchedIds::const_iterator iRecH = matchedIds.begin(); iRecH != matchedIds.end(); ++iRecH) {
1375     // vector of associated simhits associated to the current rechit
1376     std::vector<SimHitIdpr> const &SimTrackIds = (*iRecH)->second;
1377 
1378     bool found = false;
1379 
1380     for (const auto &iSimH : SimTrackIds) {
1381       uint32_t simtrackId = iSimH.first;
1382       EncodedEventId evtId = iSimH.second;
1383 
1384       // look for shared hits with the given TrackingParticle (looping over
1385       // component SimTracks)
1386       for (const auto &simtrack : g4Tracks) {
1387         if (simtrack.trackId() == simtrackId && simtrack.eventId() == evtId) {
1388           found = true;
1389           break;
1390         }
1391       }
1392 
1393       if (found) {
1394         nshared++;
1395         break;
1396       }
1397     }
1398   }
1399 
1400   return nshared;
1401 }
1402 
1403 std::string MuonAssociatorByHitsHelper::write_matched_simtracks(const std::vector<SimHitIdpr> &SimTrackIds) const {
1404   if (SimTrackIds.empty())
1405     return "  *** UNMATCHED ***";
1406 
1407   string hitlog(" matched to SimTrack");
1408 
1409   for (size_t j = 0; j < SimTrackIds.size(); j++) {
1410     char buf[64];
1411     snprintf(buf,
1412              64,
1413              " Id:%i/Evt:(%i,%i) ",
1414              SimTrackIds[j].first,
1415              SimTrackIds[j].second.event(),
1416              SimTrackIds[j].second.bunchCrossing());
1417     hitlog += buf;
1418   }
1419   return hitlog;
1420 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!