Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-05-10 05:44:19

0001 #include <iostream>
0002 #include <fstream>
0003 #include <type_traits>
0004 #include <typeinfo>
0005 
0006 #include "QuickTrackAssociatorByHitsImpl.h"
0007 
0008 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0009 
0010 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "DataFormats/Math/interface/deltaR.h"
0013 
0014 #include "SimTracker/TrackAssociation/interface/trackHitsToClusterRefs.h"
0015 
0016 //
0017 // Use the unnamed namespace for utility functions only used in this file
0018 //
0019 namespace {
0020   //
0021   // All of these functions are pretty straightforward but the implementation is type dependent.
0022   // The templated methods call these for the type specific parts and the compiler will resolve
0023   // the type and call the correct overload.
0024   //
0025 
0026   template <class T_element>
0027   size_t collectionSize(const edm::RefToBaseVector<T_element>& collection) {
0028     return collection.size();
0029   }
0030 
0031   template <class T_element>
0032   size_t collectionSize(const edm::Handle<T_element>& hCollection) {
0033     return hCollection->size();
0034   }
0035 
0036   template <class T_element>
0037   size_t collectionSize(const edm::RefVector<T_element>& collection) {
0038     return collection.size();
0039   }
0040 
0041   template <typename T_Track>
0042   const T_Track* getTrackAt(const edm::RefToBaseVector<T_Track>& trackCollection, size_t index) {
0043     return &*trackCollection[index];  // pretty obscure dereference
0044   }
0045 
0046   template <typename T_Coll>
0047   const auto* getTrackAt(const edm::Handle<T_Coll>& pTrackCollection, size_t index) {
0048     return &(*pTrackCollection)[index];
0049   }
0050 
0051   template <typename T_Track>
0052   int numberOfValidHits(const T_Track& track) {
0053     return std::count_if(track.recHits().begin(), track.recHits().end(), [](auto const& h) { return h.isValid(); });
0054   }
0055   template <>
0056   int numberOfValidHits(const reco::Track& track) {
0057     return track.numberOfValidHits();
0058   }
0059 
0060   const TrackingParticle* getTrackingParticleAt(const edm::Handle<TrackingParticleCollection>& pCollection,
0061                                                 size_t index) {
0062     return &(*pCollection.product())[index];
0063   }
0064 
0065   const TrackingParticle* getTrackingParticleAt(const edm::RefVector<TrackingParticleCollection>& collection,
0066                                                 size_t index) {
0067     return &*collection[index];
0068   }
0069 
0070   template <typename T_Track>
0071   edm::RefToBase<T_Track> getRefToTrackAt(const edm::RefToBaseVector<T_Track>& trackCollection, size_t index) {
0072     return trackCollection[index];
0073   }
0074   template <typename T_Track>
0075   edm::RefToBase<T_Track> getRefToTrackAt(const edm::Handle<edm::View<T_Track>>& pTrackCollection, size_t index) {
0076     return edm::RefToBase<T_Track>(pTrackCollection, index);
0077   }
0078   template <typename T_Track>
0079   edm::Ref<std::vector<T_Track>> getRefToTrackAt(const edm::Handle<std::vector<T_Track>>& pTrackCollection,
0080                                                  size_t index) {
0081     return edm::Ref<std::vector<T_Track>>(pTrackCollection, index);
0082   }
0083 
0084   edm::Ref<TrackingParticleCollection> getRefToTrackingParticleAt(
0085       const edm::Handle<TrackingParticleCollection>& pCollection, size_t index) {
0086     return edm::Ref<TrackingParticleCollection>(pCollection, index);
0087   }
0088 
0089   edm::Ref<TrackingParticleCollection> getRefToTrackingParticleAt(
0090       const edm::RefVector<TrackingParticleCollection>& collection, size_t index) {
0091     return collection[index];
0092   }
0093 
0094   void fillKeys(edm::IndexSet& keys, const edm::RefVector<TrackingParticleCollection>& collection) {
0095     keys.reserve(collection.size());
0096     for (const auto& ref : collection) {
0097       keys.insert(ref.key());
0098     }
0099   }
0100 
0101   template <typename Coll>
0102   void checkClusterMapProductID(const ClusterTPAssociation& clusterToTPMap, const Coll& collection) {
0103     clusterToTPMap.checkMappedProductID(collection.id());
0104   }
0105 
0106   template <typename Coll>
0107   void checkClusterMapProductID(const TrackerHitAssociator& hitAssociator, const Coll& collection) {}
0108 
0109 }  // namespace
0110 
0111 QuickTrackAssociatorByHitsImpl::QuickTrackAssociatorByHitsImpl(edm::EDProductGetter const& productGetter,
0112                                                                std::unique_ptr<const TrackerHitAssociator> hitAssoc,
0113                                                                const ClusterTPAssociation* clusterToTPMap,
0114 
0115                                                                bool absoluteNumberOfHits,
0116                                                                double qualitySimToReco,
0117                                                                double puritySimToReco,
0118                                                                double pixelHitWeight,
0119                                                                double cutRecoToSim,
0120                                                                bool threeHitTracksAreSpecial,
0121                                                                SimToRecoDenomType simToRecoDenominator)
0122     : productGetter_(&productGetter),
0123       hitAssociator_(std::move(hitAssoc)),
0124       clusterToTPMap_(clusterToTPMap),
0125       qualitySimToReco_(qualitySimToReco),
0126       puritySimToReco_(puritySimToReco),
0127       pixelHitWeight_(pixelHitWeight),
0128       cutRecoToSim_(cutRecoToSim),
0129       simToRecoDenominator_(simToRecoDenominator),
0130       threeHitTracksAreSpecial_(threeHitTracksAreSpecial),
0131       absoluteNumberOfHits_(absoluteNumberOfHits) {}
0132 
0133 reco::RecoToSimCollection QuickTrackAssociatorByHitsImpl::associateRecoToSim(
0134     const edm::Handle<edm::View<reco::Track>>& trackCollectionHandle,
0135     const edm::Handle<TrackingParticleCollection>& trackingParticleCollectionHandle) const {
0136   // Only pass the one that was successfully created to the templated method.
0137   if (not clusterToTPMap_)
0138     return associateRecoToSimImplementation<reco::RecoToSimCollection>(
0139         trackCollectionHandle, trackingParticleCollectionHandle, nullptr, *hitAssociator_);
0140   else
0141     return associateRecoToSimImplementation<reco::RecoToSimCollection>(
0142         trackCollectionHandle, trackingParticleCollectionHandle, nullptr, *clusterToTPMap_);
0143 }
0144 
0145 reco::SimToRecoCollection QuickTrackAssociatorByHitsImpl::associateSimToReco(
0146     const edm::Handle<edm::View<reco::Track>>& trackCollectionHandle,
0147     const edm::Handle<TrackingParticleCollection>& trackingParticleCollectionHandle) const {
0148   // Only pass the one that was successfully created to the templated method.
0149   if (not clusterToTPMap_)
0150     return associateSimToRecoImplementation<reco::SimToRecoCollection>(
0151         trackCollectionHandle, trackingParticleCollectionHandle, nullptr, *hitAssociator_);
0152   else
0153     return associateSimToRecoImplementation<reco::SimToRecoCollection>(
0154         trackCollectionHandle, trackingParticleCollectionHandle, nullptr, *clusterToTPMap_);
0155 }
0156 
0157 reco::RecoToSimCollection QuickTrackAssociatorByHitsImpl::associateRecoToSim(
0158     const edm::RefToBaseVector<reco::Track>& trackCollection,
0159     const edm::RefVector<TrackingParticleCollection>& trackingParticleCollection) const {
0160   // Only pass the one that was successfully created to the templated method.
0161   if (not clusterToTPMap_)
0162     return associateRecoToSimImplementation<reco::RecoToSimCollection>(
0163         trackCollection, trackingParticleCollection, nullptr, *hitAssociator_);
0164   else {
0165     TrackingParticleRefKeySet tpKeys;
0166     fillKeys(tpKeys, trackingParticleCollection);
0167     return associateRecoToSimImplementation<reco::RecoToSimCollection>(
0168         trackCollection, trackingParticleCollection, &tpKeys, *clusterToTPMap_);
0169   }
0170 }
0171 
0172 reco::SimToRecoCollection QuickTrackAssociatorByHitsImpl::associateSimToReco(
0173     const edm::RefToBaseVector<reco::Track>& trackCollection,
0174     const edm::RefVector<TrackingParticleCollection>& trackingParticleCollection) const {
0175   // Only pass the one that was successfully created to the templated method.
0176   if (not clusterToTPMap_)
0177     return associateSimToRecoImplementation<reco::SimToRecoCollection>(
0178         trackCollection, trackingParticleCollection, nullptr, *hitAssociator_);
0179   else {
0180     TrackingParticleRefKeySet tpKeys;
0181     fillKeys(tpKeys, trackingParticleCollection);
0182     return associateSimToRecoImplementation<reco::SimToRecoCollection>(
0183         trackCollection, trackingParticleCollection, &tpKeys, *clusterToTPMap_);
0184   }
0185 }
0186 
0187 template <class T_RecoToSimCollection,
0188           class T_TrackCollection,
0189           class T_TrackingParticleCollection,
0190           class T_hitOrClusterAssociator>
0191 T_RecoToSimCollection QuickTrackAssociatorByHitsImpl::associateRecoToSimImplementation(
0192     const T_TrackCollection& trackCollection,
0193     const T_TrackingParticleCollection& trackingParticleCollection,
0194     const TrackingParticleRefKeySet* trackingParticleKeys,
0195     T_hitOrClusterAssociator hitOrClusterAssociator) const {
0196   T_RecoToSimCollection returnValue(productGetter_);
0197   if (::collectionSize(trackingParticleCollection) == 0)
0198     return returnValue;
0199 
0200   checkClusterMapProductID(hitOrClusterAssociator, trackingParticleCollection);
0201 
0202   size_t collectionSize = ::collectionSize(trackCollection);  // Delegate away type specific part
0203 
0204   for (size_t i = 0; i < collectionSize; ++i) {
0205     // Get a normal pointer for ease of use. This part is type specific so delegate.
0206     const auto* pTrack = ::getTrackAt(trackCollection, i);
0207 
0208     // The return of this function has first as the index and second as the number of associated hits
0209     std::vector<std::pair<edm::Ref<TrackingParticleCollection>, double>> trackingParticleQualityPairs =
0210         associateTrack(hitOrClusterAssociator,
0211                        trackingParticleCollection,
0212                        trackingParticleKeys,
0213                        pTrack->recHits().begin(),
0214                        pTrack->recHits().end());
0215 
0216     // int nt = 0;
0217     for (auto iTrackingParticleQualityPair = trackingParticleQualityPairs.begin();
0218          iTrackingParticleQualityPair != trackingParticleQualityPairs.end();
0219          ++iTrackingParticleQualityPair) {
0220       const edm::Ref<TrackingParticleCollection>& trackingParticleRef = iTrackingParticleQualityPair->first;
0221       double numberOfSharedClusters = iTrackingParticleQualityPair->second;
0222       double numberOfValidTrackClusters = weightedNumberOfTrackClusters(*pTrack, hitOrClusterAssociator);
0223 
0224       if (numberOfSharedClusters == 0.0)
0225         continue;  // No point in continuing if there was no association
0226 
0227       //if electron subtract double counting
0228       if (abs(trackingParticleRef->pdgId()) == 11 &&
0229           (trackingParticleRef->g4Track_end() - trackingParticleRef->g4Track_begin()) > 1) {
0230         numberOfSharedClusters -= getDoubleCount(
0231             hitOrClusterAssociator, pTrack->recHits().begin(), pTrack->recHits().end(), trackingParticleRef);
0232       }
0233 
0234       double quality;
0235       if (absoluteNumberOfHits_)
0236         quality = numberOfSharedClusters;
0237       else if (numberOfValidTrackClusters != 0.0)
0238         quality = numberOfSharedClusters / numberOfValidTrackClusters;
0239       else
0240         quality = 0;
0241       if (quality > cutRecoToSim_ &&
0242           !(threeHitTracksAreSpecial_ && ::numberOfValidHits(*pTrack) == 3 && numberOfSharedClusters < 3.0)) {
0243         // Getting the RefToBase is dependent on the type of trackCollection, so delegate that to an overload.
0244         returnValue.insert(::getRefToTrackAt(trackCollection, i), std::make_pair(trackingParticleRef, quality));
0245       }
0246     }
0247   }
0248   returnValue.post_insert();
0249   return returnValue;
0250 }
0251 
0252 template <class T_SimToRecoCollection,
0253           class T_TrackCollection,
0254           class T_TrackingParticleCollection,
0255           class T_hitOrClusterAssociator>
0256 T_SimToRecoCollection QuickTrackAssociatorByHitsImpl::associateSimToRecoImplementation(
0257     const T_TrackCollection& trackCollection,
0258     const T_TrackingParticleCollection& trackingParticleCollection,
0259     const TrackingParticleRefKeySet* trackingParticleKeys,
0260     T_hitOrClusterAssociator hitOrClusterAssociator) const {
0261   T_SimToRecoCollection returnValue(productGetter_);
0262   if (::collectionSize(trackingParticleCollection) == 0)
0263     return returnValue;
0264 
0265   checkClusterMapProductID(hitOrClusterAssociator, trackingParticleCollection);
0266 
0267   size_t collectionSize = ::collectionSize(trackCollection);  // Delegate away type specific part
0268 
0269   for (size_t i = 0; i < collectionSize; ++i) {
0270     // Get a normal pointer for ease of use. This part is type specific so delegate.
0271     const auto* pTrack = ::getTrackAt(trackCollection, i);
0272 
0273     // The return of this function has first as an edm:Ref to the associated TrackingParticle, and second as the number of associated hits
0274     std::vector<std::pair<edm::Ref<TrackingParticleCollection>, double>> trackingParticleQualityPairs =
0275         associateTrack(hitOrClusterAssociator,
0276                        trackingParticleCollection,
0277                        trackingParticleKeys,
0278                        pTrack->recHits().begin(),
0279                        pTrack->recHits().end());
0280 
0281     // int nt = 0;
0282     for (auto iTrackingParticleQualityPair = trackingParticleQualityPairs.begin();
0283          iTrackingParticleQualityPair != trackingParticleQualityPairs.end();
0284          ++iTrackingParticleQualityPair) {
0285       const edm::Ref<TrackingParticleCollection>& trackingParticleRef = iTrackingParticleQualityPair->first;
0286       double numberOfSharedClusters = iTrackingParticleQualityPair->second;
0287       double numberOfValidTrackClusters = weightedNumberOfTrackClusters(*pTrack, hitOrClusterAssociator);
0288       size_t numberOfSimulatedHits = 0;  // Set a few lines below, but only if required.
0289 
0290       if (numberOfSharedClusters == 0.0)
0291         continue;  // No point in continuing if there was no association
0292 
0293       if (simToRecoDenominator_ == denomsim ||
0294           (numberOfSharedClusters < 3.0 &&
0295            threeHitTracksAreSpecial_))  // the numberOfSimulatedHits is not always required, so can skip counting in some circumstances
0296       {
0297         // Note that in the standard TrackAssociatorByHits, all of the hits in associatedTrackingParticleHits are checked for
0298         // various things.  I'm not sure what these checks are for but they depend on the UseGrouping and UseSplitting settings.
0299         // This associator works as though both UseGrouping and UseSplitting were set to true, i.e. just counts the number of
0300         // hits in the tracker.
0301         numberOfSimulatedHits = trackingParticleRef->numberOfTrackerHits();
0302       }
0303 
0304       //if electron subtract double counting
0305       if (abs(trackingParticleRef->pdgId()) == 11 &&
0306           (trackingParticleRef->g4Track_end() - trackingParticleRef->g4Track_begin()) > 1) {
0307         numberOfSharedClusters -= getDoubleCount(
0308             hitOrClusterAssociator, pTrack->recHits().begin(), pTrack->recHits().end(), trackingParticleRef);
0309       }
0310 
0311       double purity = numberOfSharedClusters / numberOfValidTrackClusters;
0312       double quality;
0313       if (absoluteNumberOfHits_)
0314         quality = numberOfSharedClusters;
0315       else if (simToRecoDenominator_ == denomsim && numberOfSimulatedHits != 0)
0316         quality = numberOfSharedClusters / static_cast<double>(numberOfSimulatedHits);
0317       else if (simToRecoDenominator_ == denomreco && numberOfValidTrackClusters != 0)
0318         quality = purity;
0319       else
0320         quality = 0;
0321 
0322       if (quality > qualitySimToReco_ &&
0323           !(threeHitTracksAreSpecial_ && numberOfSimulatedHits == 3 && numberOfSharedClusters < 3.0) &&
0324           (absoluteNumberOfHits_ || (purity > puritySimToReco_))) {
0325         // Getting the RefToBase is dependent on the type of trackCollection, so delegate that to an overload.
0326         returnValue.insert(trackingParticleRef, std::make_pair(::getRefToTrackAt(trackCollection, i), quality));
0327       }
0328     }
0329   }
0330   returnValue.post_insert();
0331   return returnValue;
0332 }
0333 
0334 template <typename T_TPCollection, typename iter>
0335 std::vector<std::pair<edm::Ref<TrackingParticleCollection>, double>> QuickTrackAssociatorByHitsImpl::associateTrack(
0336     const TrackerHitAssociator& hitAssociator,
0337     const T_TPCollection& trackingParticles,
0338     const TrackingParticleRefKeySet* trackingParticleKeys,
0339     iter begin,
0340     iter end) const {
0341   // The pairs in this vector have a Ref to the associated TrackingParticle as "first" and the weighted number of associated hits as "second"
0342   std::vector<std::pair<edm::Ref<TrackingParticleCollection>, double>> returnValue;
0343 
0344   // The pairs in this vector have first as the sim track identifiers, and second the number of reco hits associated to that sim track.
0345   // Most reco hits will probably have come from the same sim track, so the number of entries in this vector should be fewer than the
0346   // number of reco hits.  The pair::second entries should add up to the total number of reco hits though.
0347   std::vector<std::pair<SimTrackIdentifiers, double>> hitIdentifiers =
0348       getAllSimTrackIdentifiers(hitAssociator, begin, end);
0349 
0350   // Loop over the TrackingParticles
0351   size_t collectionSize = ::collectionSize(trackingParticles);
0352 
0353   for (size_t i = 0; i < collectionSize; ++i) {
0354     const TrackingParticle* pTrackingParticle = getTrackingParticleAt(trackingParticles, i);
0355 
0356     // Historically there was a requirement that pTrackingParticle->numberOfHits() > 0
0357     // However, in TrackingTruthAccumulator, the numberOfHits is calculated from a subset
0358     // of the SimHits of the SimTracks of a TrackingParticle (essentially limiting the
0359     // processType and particleType to those of the "first" hit, and particleType to the pdgId of the SimTrack).
0360     // But, here the association between tracks and TrackingParticles is done with *all* the hits of
0361     // TrackingParticle, so we should not rely on the numberOfHits() calculated with a subset of SimHits.
0362 
0363     double numberOfAssociatedHits = 0;
0364     // Loop over all of the sim track identifiers and see if any of them are part of this TrackingParticle. If they are, add
0365     // the number of reco hits associated to that sim track to the total number of associated hits.
0366     for (const auto& identifierCountPair : hitIdentifiers) {
0367       if (trackingParticleContainsIdentifier(pTrackingParticle, identifierCountPair.first))
0368         numberOfAssociatedHits += identifierCountPair.second;
0369     }
0370 
0371     if (numberOfAssociatedHits > 0) {
0372       returnValue.push_back(std::make_pair(getRefToTrackingParticleAt(trackingParticles, i), numberOfAssociatedHits));
0373     }
0374   }
0375 
0376   return returnValue;
0377 }
0378 
0379 template <typename T_TPCollection, typename iter>
0380 std::vector<std::pair<edm::Ref<TrackingParticleCollection>, double>> QuickTrackAssociatorByHitsImpl::associateTrack(
0381     const ClusterTPAssociation& clusterToTPMap,
0382     const T_TPCollection& trackingParticles,
0383     const TrackingParticleRefKeySet* trackingParticleKeys,
0384     iter begin,
0385     iter end) const {
0386   // Note that the trackingParticles parameter is not actually required since all the information is in clusterToTPMap,
0387   // but the method signature has to match the other overload because it is called from a templated method.
0388 
0389   // Note further, that we can't completely ignore the
0390   // trackingParticles parameter, in case it is a subset of those
0391   // TrackingParticles used to construct clusterToTPMap (via the
0392   // TrackingParticleRefVector overloads). The trackingParticles
0393   // parameter is still ignored since looping over it on every call
0394   // would be expensive, but the keys of the TrackingParticleRefs are
0395   // cached to an IndexSet (trackingParticleKeys) which is used
0396   // as a fast search structure.
0397 
0398   // The pairs in this vector have a Ref to the associated TrackingParticle as "first" and the weighted number of associated clusters as "second"
0399   // Note: typedef edm::Ref<TrackingParticleCollection> TrackingParticleRef;
0400   std::vector<std::pair<edm::Ref<TrackingParticleCollection>, double>> returnValue;
0401   if (clusterToTPMap.empty())
0402     return returnValue;
0403 
0404   // The pairs in this vector have first as the TP, and second the number of reco clusters associated to that TP.
0405   // Most reco clusters will probably have come from the same sim track (i.e TP), so the number of entries in this
0406   // vector should be fewer than the number of clusters. The pair::second entries should add up to the total
0407   // number of reco clusters though.
0408   std::vector<OmniClusterRef> oClusters = track_associator::hitsToClusterRefs(begin, end);
0409 
0410   std::map<TrackingParticleRef, double> lmap;
0411   for (std::vector<OmniClusterRef>::const_iterator it = oClusters.begin(); it != oClusters.end(); ++it) {
0412     auto range = clusterToTPMap.equal_range(*it);
0413     const double weight = it->isPixel() ? pixelHitWeight_ : 1.0;
0414     if (range.first != range.second) {
0415       for (auto ip = range.first; ip != range.second; ++ip) {
0416         const TrackingParticleRef trackingParticle = (ip->second);
0417 
0418         if (trackingParticleKeys && !trackingParticleKeys->has(trackingParticle.key()))
0419           continue;
0420 
0421         // Historically there was a requirement that pTrackingParticle->numberOfHits() > 0
0422         // However, in TrackingTruthAccumulator, the numberOfHits is calculated from a subset
0423         // of the SimHits of the SimTracks of a TrackingParticle (essentially limiting the
0424         // processType and particleType to those of the "first" hit, and particleType to the pdgId of the SimTrack).
0425         // But, here the association between tracks and TrackingParticles is done with *all* the hits of
0426         // TrackingParticle, so we should not rely on the numberOfHits() calculated with a subset of SimHits.
0427 
0428         /* Alternative implementation to avoid the use of lmap... memory slightly improved but slightly slower...
0429                  std::pair<edm::Ref<TrackingParticleCollection>,size_t> tpIntPair(trackingParticle, 1);
0430                  auto tp_range = std::equal_range(returnValue.begin(), returnValue.end(), tpIntPair, tpIntPairGreater);
0431                  if ((tp_range.second-tp_range.first)>1) {
0432                  edm::LogError("TrackAssociator") << ">>> Error in counting TPs!" << " file: " << __FILE__ << " line: " << __LINE__;
0433                  }
0434                  if(tp_range.first != tp_range.second) {
0435                  tp_range.first->second++;
0436                  } else {
0437                  returnValue.push_back(tpIntPair);
0438                  std::sort(returnValue.begin(), returnValue.end(), tpIntPairGreater);
0439                  }
0440                  */
0441         auto jpos = lmap.find(trackingParticle);
0442         if (jpos != lmap.end())
0443           jpos->second += weight;
0444         else
0445           lmap.insert(std::make_pair(trackingParticle, weight));
0446       }
0447     }
0448   }
0449   // now copy the map to returnValue
0450   for (auto ip = lmap.begin(); ip != lmap.end(); ++ip) {
0451     returnValue.push_back(std::make_pair(ip->first, ip->second));
0452   }
0453   return returnValue;
0454 }
0455 
0456 template <typename iter>
0457 std::vector<std::pair<QuickTrackAssociatorByHitsImpl::SimTrackIdentifiers, double>>
0458 QuickTrackAssociatorByHitsImpl::getAllSimTrackIdentifiers(const TrackerHitAssociator& hitAssociator,
0459                                                           iter begin,
0460                                                           iter end) const {
0461   // The pairs in this vector have first as the sim track identifiers, and second the number of reco hits associated to that sim track.
0462   std::vector<std::pair<SimTrackIdentifiers, double>> returnValue;
0463 
0464   std::vector<SimTrackIdentifiers> simTrackIdentifiers;
0465   // Loop over all of the rec hits in the track
0466   //iter tRHIterBeginEnd = getTRHIterBeginEnd( pTrack );
0467   for (iter iRecHit = begin; iRecHit != end; ++iRecHit) {
0468     if (track_associator::getHitFromIter(iRecHit)->isValid()) {
0469       simTrackIdentifiers.clear();
0470 
0471       // Get the identifiers for the sim track that this hit came from. There should only be one entry unless clusters
0472       // have merged (as far as I know).
0473       hitAssociator.associateHitId(*(track_associator::getHitFromIter(iRecHit)),
0474                                    simTrackIdentifiers);  // This call fills simTrackIdentifiers
0475 
0476       const auto subdetId = track_associator::getHitFromIter(iRecHit)->geographicalId().subdetId();
0477       const double weight = (subdetId == PixelSubdetector::PixelBarrel || subdetId == PixelSubdetector::PixelEndcap)
0478                                 ? pixelHitWeight_
0479                                 : 1.0;
0480 
0481       // Loop over each identifier, and add it to the return value only if it's not already in there
0482       for (std::vector<SimTrackIdentifiers>::const_iterator iIdentifier = simTrackIdentifiers.begin();
0483            iIdentifier != simTrackIdentifiers.end();
0484            ++iIdentifier) {
0485         std::vector<std::pair<SimTrackIdentifiers, double>>::iterator iIdentifierCountPair;
0486         for (auto iIdentifierCountPair = returnValue.begin(); iIdentifierCountPair != returnValue.end();
0487              ++iIdentifierCountPair) {
0488           if (iIdentifierCountPair->first.first == iIdentifier->first &&
0489               iIdentifierCountPair->first.second == iIdentifier->second) {
0490             // This sim track identifier is already in the list, so increment the count of how many hits it relates to.
0491             iIdentifierCountPair->second += weight;
0492             break;
0493           }
0494         }
0495         if (iIdentifierCountPair == returnValue.end())
0496           returnValue.push_back(std::make_pair(*iIdentifier, 1.0));
0497         // This identifier wasn't found, so add it
0498       }
0499     }
0500   }
0501   return returnValue;
0502 }
0503 
0504 bool QuickTrackAssociatorByHitsImpl::trackingParticleContainsIdentifier(const TrackingParticle* pTrackingParticle,
0505                                                                         const SimTrackIdentifiers& identifier) const {
0506   // Loop over all of the g4 tracks in the tracking particle
0507   for (std::vector<SimTrack>::const_iterator iSimTrack = pTrackingParticle->g4Track_begin();
0508        iSimTrack != pTrackingParticle->g4Track_end();
0509        ++iSimTrack) {
0510     // And see if the sim track identifiers match
0511     if (iSimTrack->eventId() == identifier.second && iSimTrack->trackId() == identifier.first) {
0512       return true;
0513     }
0514   }
0515 
0516   // If control has made it this far then none of the identifiers were found in
0517   // any of the g4 tracks, so return false.
0518   return false;
0519 }
0520 
0521 template <typename iter>
0522 double QuickTrackAssociatorByHitsImpl::getDoubleCount(const TrackerHitAssociator& hitAssociator,
0523                                                       iter startIterator,
0524                                                       iter endIterator,
0525                                                       TrackingParticleRef associatedTrackingParticle) const {
0526   // This method is largely copied from the standard TrackAssociatorByHits. Once I've tested how much difference
0527   // it makes I'll go through and comment it properly.
0528 
0529   // FIXME: It may be that this piece is not fully correct for
0530   // counting how many times a single *cluster* is matched to many
0531   // SimTracks of a single TrackingParticle (see comments in
0532   // getDoubleCount(ClusterTPAssociation) overload). To be verified
0533   // some time.
0534 
0535   double doubleCount = 0.0;
0536   std::vector<SimHitIdpr> SimTrackIdsDC;
0537 
0538   for (iter iHit = startIterator; iHit != endIterator; iHit++) {
0539     int idcount = 0;
0540 
0541     SimTrackIdsDC.clear();
0542     hitAssociator.associateHitId(*(track_associator::getHitFromIter(iHit)), SimTrackIdsDC);
0543     if (SimTrackIdsDC.size() > 1) {
0544       for (TrackingParticle::g4t_iterator g4T = associatedTrackingParticle->g4Track_begin();
0545            g4T != associatedTrackingParticle->g4Track_end();
0546            ++g4T) {
0547         if (find(SimTrackIdsDC.begin(),
0548                  SimTrackIdsDC.end(),
0549                  SimHitIdpr((*g4T).trackId(), SimTrackIdsDC.begin()->second)) != SimTrackIdsDC.end()) {
0550           idcount++;
0551         }
0552       }
0553     }
0554     if (idcount > 1) {
0555       const auto subdetId = track_associator::getHitFromIter(iHit)->geographicalId().subdetId();
0556       const double weight = (subdetId == PixelSubdetector::PixelBarrel || subdetId == PixelSubdetector::PixelEndcap)
0557                                 ? pixelHitWeight_
0558                                 : 1.0;
0559       doubleCount += weight * (idcount - 1);
0560     }
0561   }
0562 
0563   return doubleCount;
0564 }
0565 
0566 template <typename iter>
0567 double QuickTrackAssociatorByHitsImpl::getDoubleCount(const ClusterTPAssociation& clusterToTPList,
0568                                                       iter startIterator,
0569                                                       iter endIterator,
0570                                                       TrackingParticleRef associatedTrackingParticle) const {
0571   // This code here was written by Subir Sarkar. I'm just splitting it off into a
0572   // separate method. - Grimes 01/May/2014
0573 
0574   // The point here is that the electron TrackingParticles may contain
0575   // multiple SimTracks (from the bremsstrahling), and (historically)
0576   // the each matched hit/cluster has been multiplied by "how many
0577   // SimTracks from the TrackingParticle" it contains charge from.
0578   // Here the amount of this double counting is calculated, so that it
0579   // can be subtracted by the calling code.
0580   //
0581   // Note that recently (hence "historically" in the paragraph above)
0582   // the ClusterTPAssociationProducer was changed to remove the
0583   // duplicate cluster->TP associations (hence making this function
0584   // obsolete), but there is more recent proof that there is some
0585   // duplication left (to be investigated).
0586 
0587   double doubleCount = 0;
0588   std::vector<SimHitIdpr> SimTrackIdsDC;
0589 
0590   for (iter iHit = startIterator; iHit != endIterator; iHit++) {
0591     std::vector<OmniClusterRef> oClusters =
0592         track_associator::hitsToClusterRefs(iHit, iHit + 1);  //only for the cluster being checked
0593     for (std::vector<OmniClusterRef>::const_iterator it = oClusters.begin(); it != oClusters.end(); ++it) {
0594       int idcount = 0;
0595 
0596       auto range = clusterToTPList.equal_range(*it);
0597       if (range.first != range.second) {
0598         for (auto ip = range.first; ip != range.second; ++ip) {
0599           const TrackingParticleRef trackingParticle = (ip->second);
0600           if (associatedTrackingParticle == trackingParticle) {
0601             idcount++;
0602           }
0603         }
0604       }
0605 
0606       if (idcount > 1) {
0607         const auto subdetId = track_associator::getHitFromIter(iHit)->geographicalId().subdetId();
0608         const double weight = (subdetId == PixelSubdetector::PixelBarrel || subdetId == PixelSubdetector::PixelEndcap)
0609                                   ? pixelHitWeight_
0610                                   : 1.0;
0611         doubleCount += weight * (idcount - 1);
0612       }
0613     }
0614   }
0615 
0616   return doubleCount;
0617 }
0618 
0619 reco::RecoToSimCollectionSeed QuickTrackAssociatorByHitsImpl::associateRecoToSim(
0620     const edm::Handle<edm::View<TrajectorySeed>>& pSeedCollectionHandle_,
0621     const edm::Handle<TrackingParticleCollection>& trackingParticleCollectionHandle) const {
0622   edm::LogVerbatim("TrackAssociator") << "Starting TrackAssociatorByHitsImpl::associateRecoToSim - #seeds="
0623                                       << pSeedCollectionHandle_->size()
0624                                       << " #TPs=" << trackingParticleCollectionHandle->size();
0625 
0626   reco::RecoToSimCollectionSeed returnValue(productGetter_);
0627 
0628   size_t collectionSize = pSeedCollectionHandle_->size();
0629 
0630   for (size_t i = 0; i < collectionSize; ++i) {
0631     const TrajectorySeed* pSeed = &(*pSeedCollectionHandle_)[i];
0632 
0633     // The return of this function has first as the index and second as the number of associated hits
0634     std::vector<std::pair<edm::Ref<TrackingParticleCollection>, double>> trackingParticleQualityPairs =
0635         (clusterToTPMap_) ? associateTrack(*clusterToTPMap_,
0636                                            trackingParticleCollectionHandle,
0637                                            nullptr,
0638                                            pSeed->recHits().begin(),
0639                                            pSeed->recHits().end())
0640                           : associateTrack(*hitAssociator_,
0641                                            trackingParticleCollectionHandle,
0642                                            nullptr,
0643                                            pSeed->recHits().begin(),
0644                                            pSeed->recHits().end());
0645     for (auto iTrackingParticleQualityPair = trackingParticleQualityPairs.begin();
0646          iTrackingParticleQualityPair != trackingParticleQualityPairs.end();
0647          ++iTrackingParticleQualityPair) {
0648       const edm::Ref<TrackingParticleCollection>& trackingParticleRef = iTrackingParticleQualityPair->first;
0649       double numberOfSharedClusters = iTrackingParticleQualityPair->second;
0650       double numberOfValidTrackClusters = clusterToTPMap_ ? weightedNumberOfTrackClusters(*pSeed, *clusterToTPMap_)
0651                                                           : weightedNumberOfTrackClusters(*pSeed, *hitAssociator_);
0652 
0653       if (numberOfSharedClusters == 0.0)
0654         continue;  // No point in continuing if there was no association
0655 
0656       //if electron subtract double counting
0657       if (abs(trackingParticleRef->pdgId()) == 11 &&
0658           (trackingParticleRef->g4Track_end() - trackingParticleRef->g4Track_begin()) > 1) {
0659         if (clusterToTPMap_)
0660           numberOfSharedClusters -=
0661               getDoubleCount(*clusterToTPMap_, pSeed->recHits().begin(), pSeed->recHits().end(), trackingParticleRef);
0662         else
0663           numberOfSharedClusters -=
0664               getDoubleCount(*hitAssociator_, pSeed->recHits().begin(), pSeed->recHits().end(), trackingParticleRef);
0665       }
0666 
0667       double quality;
0668       if (absoluteNumberOfHits_)
0669         quality = numberOfSharedClusters;
0670       else if (numberOfValidTrackClusters != 0.0)
0671         quality = numberOfSharedClusters / numberOfValidTrackClusters;
0672       else
0673         quality = 0;
0674 
0675       if (quality > cutRecoToSim_ &&
0676           !(threeHitTracksAreSpecial_ && pSeed->nHits() == 3 && numberOfSharedClusters < 3.0)) {
0677         returnValue.insert(edm::RefToBase<TrajectorySeed>(pSeedCollectionHandle_, i),
0678                            std::make_pair(trackingParticleRef, quality));
0679       }
0680     }
0681   }
0682 
0683   LogTrace("TrackAssociator") << "% of Assoc Seeds="
0684                               << ((double)returnValue.size()) / ((double)pSeedCollectionHandle_->size());
0685   returnValue.post_insert();
0686   return returnValue;
0687 }
0688 
0689 reco::SimToRecoCollectionSeed QuickTrackAssociatorByHitsImpl::associateSimToReco(
0690     const edm::Handle<edm::View<TrajectorySeed>>& pSeedCollectionHandle_,
0691     const edm::Handle<TrackingParticleCollection>& trackingParticleCollectionHandle) const {
0692   edm::LogVerbatim("TrackAssociator") << "Starting TrackAssociatorByHitsImpl::associateSimToReco - #seeds="
0693                                       << pSeedCollectionHandle_->size()
0694                                       << " #TPs=" << trackingParticleCollectionHandle->size();
0695 
0696   reco::SimToRecoCollectionSeed returnValue(productGetter_);
0697   if (trackingParticleCollectionHandle->empty())
0698     return returnValue;
0699 
0700   if (clusterToTPMap_) {
0701     checkClusterMapProductID(*clusterToTPMap_, trackingParticleCollectionHandle);
0702   }
0703 
0704   size_t collectionSize = pSeedCollectionHandle_->size();
0705 
0706   for (size_t i = 0; i < collectionSize; ++i) {
0707     const TrajectorySeed* pSeed = &(*pSeedCollectionHandle_)[i];
0708 
0709     // The return of this function has first as an edm:Ref to the associated TrackingParticle, and second as the number of associated hits
0710     std::vector<std::pair<edm::Ref<TrackingParticleCollection>, double>> trackingParticleQualityPairs =
0711         (clusterToTPMap_) ? associateTrack(*clusterToTPMap_,
0712                                            trackingParticleCollectionHandle,
0713                                            nullptr,
0714                                            pSeed->recHits().begin(),
0715                                            pSeed->recHits().end())
0716                           : associateTrack(*hitAssociator_,
0717                                            trackingParticleCollectionHandle,
0718                                            nullptr,
0719                                            pSeed->recHits().begin(),
0720                                            pSeed->recHits().end());
0721     for (auto iTrackingParticleQualityPair = trackingParticleQualityPairs.begin();
0722          iTrackingParticleQualityPair != trackingParticleQualityPairs.end();
0723          ++iTrackingParticleQualityPair) {
0724       const edm::Ref<TrackingParticleCollection>& trackingParticleRef = iTrackingParticleQualityPair->first;
0725       double numberOfSharedClusters = iTrackingParticleQualityPair->second;
0726       double numberOfValidTrackClusters = clusterToTPMap_ ? weightedNumberOfTrackClusters(*pSeed, *clusterToTPMap_)
0727                                                           : weightedNumberOfTrackClusters(*pSeed, *hitAssociator_);
0728       size_t numberOfSimulatedHits = 0;  // Set a few lines below, but only if required.
0729 
0730       if (numberOfSharedClusters == 0.0)
0731         continue;  // No point in continuing if there was no association
0732 
0733       //if electron subtract double counting
0734       if (abs(trackingParticleRef->pdgId()) == 11 &&
0735           (trackingParticleRef->g4Track_end() - trackingParticleRef->g4Track_begin()) > 1) {
0736         if (clusterToTPMap_)
0737           numberOfSharedClusters -=
0738               getDoubleCount(*clusterToTPMap_, pSeed->recHits().begin(), pSeed->recHits().end(), trackingParticleRef);
0739         else
0740           numberOfSharedClusters -=
0741               getDoubleCount(*hitAssociator_, pSeed->recHits().begin(), pSeed->recHits().end(), trackingParticleRef);
0742       }
0743 
0744       if (simToRecoDenominator_ == denomsim ||
0745           (numberOfSharedClusters < 3.0 &&
0746            threeHitTracksAreSpecial_))  // the numberOfSimulatedHits is not always required, so can skip counting in some circumstances
0747       {
0748         // Note that in the standard TrackAssociatorByHits, all of the hits in associatedTrackingParticleHits are checked for
0749         // various things.  I'm not sure what these checks are for but they depend on the UseGrouping and UseSplitting settings.
0750         // This associator works as though both UseGrouping and UseSplitting were set to true, i.e. just counts the number of
0751         // hits in the tracker.
0752         numberOfSimulatedHits = trackingParticleRef->numberOfTrackerHits();
0753       }
0754 
0755       double purity = numberOfSharedClusters / numberOfValidTrackClusters;
0756       double quality;
0757       if (absoluteNumberOfHits_)
0758         quality = numberOfSharedClusters;
0759       else if (simToRecoDenominator_ == denomsim && numberOfSimulatedHits != 0)
0760         quality = numberOfSharedClusters / static_cast<double>(numberOfSimulatedHits);
0761       else if (simToRecoDenominator_ == denomreco && numberOfValidTrackClusters != 0.0)
0762         quality = purity;
0763       else
0764         quality = 0;
0765 
0766       if (quality > qualitySimToReco_ &&
0767           !(threeHitTracksAreSpecial_ && numberOfSimulatedHits == 3 && numberOfSharedClusters < 3.0) &&
0768           (absoluteNumberOfHits_ || (purity > puritySimToReco_))) {
0769         returnValue.insert(trackingParticleRef,
0770                            std::make_pair(edm::RefToBase<TrajectorySeed>(pSeedCollectionHandle_, i), quality));
0771       }
0772     }
0773   }
0774 
0775   LogTrace("TrackAssociator") << "% of Assoc TPs="
0776                               << ((double)returnValue.size()) / ((double)trackingParticleCollectionHandle->size());
0777   returnValue.post_insert();
0778   return returnValue;
0779 }
0780 
0781 reco::RecoToSimCollectionTCandidate QuickTrackAssociatorByHitsImpl::associateRecoToSim(
0782     const edm::Handle<TrackCandidateCollection>& trackCollectionHandle,
0783     const edm::RefVector<TrackingParticleCollection>& trackingParticleCollectionHandle) const {
0784   // Only pass the one that was successfully created to the templated method.
0785   if (not clusterToTPMap_)
0786     return associateRecoToSimImplementation<reco::RecoToSimCollectionTCandidate>(
0787         trackCollectionHandle, trackingParticleCollectionHandle, nullptr, *hitAssociator_);
0788   else
0789     return associateRecoToSimImplementation<reco::RecoToSimCollectionTCandidate>(
0790         trackCollectionHandle, trackingParticleCollectionHandle, nullptr, *clusterToTPMap_);
0791 }
0792 
0793 reco::SimToRecoCollectionTCandidate QuickTrackAssociatorByHitsImpl::associateSimToReco(
0794     const edm::Handle<TrackCandidateCollection>& trackCollectionHandle,
0795     const edm::RefVector<TrackingParticleCollection>& trackingParticleCollectionHandle) const {
0796   // Only pass the one that was successfully created to the templated method.
0797   if (not clusterToTPMap_)
0798     return associateSimToRecoImplementation<reco::SimToRecoCollectionTCandidate>(
0799         trackCollectionHandle, trackingParticleCollectionHandle, nullptr, *hitAssociator_);
0800   else
0801     return associateSimToRecoImplementation<reco::SimToRecoCollectionTCandidate>(
0802         trackCollectionHandle, trackingParticleCollectionHandle, nullptr, *clusterToTPMap_);
0803 }
0804 
0805 // count hits
0806 template <>
0807 double QuickTrackAssociatorByHitsImpl::weightedNumberOfTrackClusters(const reco::Track& track,
0808                                                                      const TrackerHitAssociator&) const {
0809   const reco::HitPattern& p = track.hitPattern();
0810   const auto pixelHits = p.numberOfValidPixelHits();
0811   const auto otherHits = p.numberOfValidHits() - pixelHits;
0812   return pixelHits * pixelHitWeight_ + otherHits;
0813 }
0814 
0815 template <typename T_Track>
0816 double QuickTrackAssociatorByHitsImpl::weightedNumberOfTrackClusters(const T_Track& seed,
0817                                                                      const TrackerHitAssociator&) const {
0818   double sum = 0.0;
0819   for (auto iHit = seed.recHits().begin(); iHit != seed.recHits().end(); ++iHit) {
0820     const auto subdetId = track_associator::getHitFromIter(iHit)->geographicalId().subdetId();
0821     const double weight = (subdetId == PixelSubdetector::PixelBarrel || subdetId == PixelSubdetector::PixelEndcap)
0822                               ? pixelHitWeight_
0823                               : 1.0;
0824     sum += weight;
0825   }
0826   return sum;
0827 }
0828 
0829 // count clusters
0830 template <>
0831 double QuickTrackAssociatorByHitsImpl::weightedNumberOfTrackClusters(const reco::Track& track,
0832                                                                      const ClusterTPAssociation&) const {
0833   return weightedNumberOfTrackClusters(track.recHitsBegin(), track.recHitsEnd());
0834 }
0835 template <typename T_Track>
0836 double QuickTrackAssociatorByHitsImpl::weightedNumberOfTrackClusters(const T_Track& seed,
0837                                                                      const ClusterTPAssociation&) const {
0838   const auto& hitRange = seed.recHits();
0839   return weightedNumberOfTrackClusters(hitRange.begin(), hitRange.end());
0840 }
0841 
0842 template <typename iter>
0843 double QuickTrackAssociatorByHitsImpl::weightedNumberOfTrackClusters(iter begin, iter end) const {
0844   double weightedClusters = 0.0;
0845   for (iter iRecHit = begin; iRecHit != end; ++iRecHit) {
0846     const auto subdetId = track_associator::getHitFromIter(iRecHit)->geographicalId().subdetId();
0847     const double weight = (subdetId == PixelSubdetector::PixelBarrel || subdetId == PixelSubdetector::PixelEndcap)
0848                               ? pixelHitWeight_
0849                               : 1.0;
0850     LogTrace("QuickTrackAssociatorByHitsImpl")
0851         << "  detId: " << track_associator::getHitFromIter(iRecHit)->geographicalId().rawId();
0852     LogTrace("QuickTrackAssociatorByHitsImpl") << "  weight: " << weight;
0853     std::vector<OmniClusterRef> oClusters =
0854         track_associator::hitsToClusterRefs(iRecHit, iRecHit + 1);  //only for the cluster being checked
0855     for (std::vector<OmniClusterRef>::const_iterator it = oClusters.begin(); it != oClusters.end(); ++it) {
0856       weightedClusters += weight;
0857     }
0858   }
0859   LogTrace("QuickTrackAssociatorByHitsImpl") << "  total weighted clusters: " << weightedClusters;
0860   return weightedClusters;
0861 }