Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix QuickTrackAssociatorByHits for composite hits #20364

Merged
merged 5 commits into from
Sep 5, 2017
Merged

Fix QuickTrackAssociatorByHits for composite hits #20364

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
287c26a
using numberOfValidTrackClusters instead of numberOfValidTrackHits
ebrondol Aug 12, 2017
963631f
consistent nomenclature
ebrondol Aug 12, 2017
1a8c051
In getDoubleCount, count the number of duplicate cluster->TP matches …
makortel Aug 16, 2017
d3e1f4a
Code style checks
makortel Sep 4, 2017
330377f
For TrackerHitAssociator count denominator as hits
makortel Sep 4, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
141 changes: 95 additions & 46 deletions SimTracker/TrackAssociatorProducers/plugins/QuickTrackAssociatorByHitsImpl.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -190,22 +190,22 @@ reco::RecoToSimCollection QuickTrackAssociatorByHitsImpl::associateRecoToSimImpl
++iTrackingParticleQualityPair )
{
const edm::Ref<TrackingParticleCollection>& trackingParticleRef=iTrackingParticleQualityPair->first;
double numberOfSharedHits=iTrackingParticleQualityPair->second;
double numberOfValidTrackHits=weightedNumberOfTrackHits(*pTrack);
double numberOfSharedClusters=iTrackingParticleQualityPair->second;
double numberOfValidTrackClusters=weightedNumberOfTrackClusters(*pTrack, hitOrClusterAssociator);

if( numberOfSharedHits == 0.0 ) continue; // No point in continuing if there was no association
if( numberOfSharedClusters == 0.0 ) continue; // No point in continuing if there was no association

//if electron subtract double counting
if( abs( trackingParticleRef->pdgId() ) == 11 && (trackingParticleRef->g4Track_end() - trackingParticleRef->g4Track_begin()) > 1 )
{
numberOfSharedHits-=getDoubleCount( hitOrClusterAssociator, pTrack->recHitsBegin(), pTrack->recHitsEnd(), trackingParticleRef );
numberOfSharedClusters-=getDoubleCount( hitOrClusterAssociator, pTrack->recHitsBegin(), pTrack->recHitsEnd(), trackingParticleRef );
}

double quality;
if( absoluteNumberOfHits_ ) quality = numberOfSharedHits;
else if( numberOfValidTrackHits != 0.0 ) quality = numberOfSharedHits / numberOfValidTrackHits;
if( absoluteNumberOfHits_ ) quality = numberOfSharedClusters;
else if( numberOfValidTrackClusters != 0.0 ) quality = numberOfSharedClusters / numberOfValidTrackClusters;
else quality=0;
if( quality > cutRecoToSim_ && !(threeHitTracksAreSpecial_ && pTrack->numberOfValidHits() == 3 && numberOfSharedHits < 3.0) )
if( quality > cutRecoToSim_ && !(threeHitTracksAreSpecial_ && pTrack->numberOfValidHits() == 3 && numberOfSharedClusters < 3.0) )
{
// Getting the RefToBase is dependent on the type of trackCollection, so delegate that to an overload.
returnValue.insert( ::getRefToTrackAt(trackCollection,i), std::make_pair( trackingParticleRef, quality ) );
Expand Down Expand Up @@ -239,13 +239,13 @@ reco::SimToRecoCollection QuickTrackAssociatorByHitsImpl::associateSimToRecoImpl
iTrackingParticleQualityPair!=trackingParticleQualityPairs.end(); ++iTrackingParticleQualityPair )
{
const edm::Ref<TrackingParticleCollection>& trackingParticleRef=iTrackingParticleQualityPair->first;
double numberOfSharedHits=iTrackingParticleQualityPair->second;
double numberOfValidTrackHits=weightedNumberOfTrackHits(*pTrack);
double numberOfSharedClusters=iTrackingParticleQualityPair->second;
double numberOfValidTrackClusters=weightedNumberOfTrackClusters(*pTrack, hitOrClusterAssociator);
size_t numberOfSimulatedHits=0; // Set a few lines below, but only if required.

if( numberOfSharedHits==0.0 ) continue; // No point in continuing if there was no association
if( numberOfSharedClusters==0.0 ) continue; // No point in continuing if there was no association

if( simToRecoDenominator_==denomsim || (numberOfSharedHits<3.0 && threeHitTracksAreSpecial_) ) // the numberOfSimulatedHits is not always required, so can skip counting in some circumstances
if( simToRecoDenominator_==denomsim || (numberOfSharedClusters<3.0 && threeHitTracksAreSpecial_) ) // the numberOfSimulatedHits is not always required, so can skip counting in some circumstances
{
// Note that in the standard TrackAssociatorByHits, all of the hits in associatedTrackingParticleHits are checked for
// various things. I'm not sure what these checks are for but they depend on the UseGrouping and UseSplitting settings.
Expand All @@ -257,17 +257,17 @@ reco::SimToRecoCollection QuickTrackAssociatorByHitsImpl::associateSimToRecoImpl
//if electron subtract double counting
if (abs(trackingParticleRef->pdgId())==11 && (trackingParticleRef->g4Track_end() - trackingParticleRef->g4Track_begin()) > 1 )
{
numberOfSharedHits -= getDoubleCount( hitOrClusterAssociator, pTrack->recHitsBegin(), pTrack->recHitsEnd(), trackingParticleRef );
numberOfSharedClusters -= getDoubleCount( hitOrClusterAssociator, pTrack->recHitsBegin(), pTrack->recHitsEnd(), trackingParticleRef );
}

double purity = numberOfSharedHits/numberOfValidTrackHits;
double purity = numberOfSharedClusters/numberOfValidTrackClusters;
double quality;
if( absoluteNumberOfHits_ ) quality = numberOfSharedHits;
else if( simToRecoDenominator_==denomsim && numberOfSimulatedHits != 0 ) quality = numberOfSharedHits/static_cast<double>(numberOfSimulatedHits);
else if( simToRecoDenominator_==denomreco && numberOfValidTrackHits != 0 ) quality=purity;
if( absoluteNumberOfHits_ ) quality = numberOfSharedClusters;
else if( simToRecoDenominator_==denomsim && numberOfSimulatedHits != 0 ) quality = numberOfSharedClusters/static_cast<double>(numberOfSimulatedHits);
else if( simToRecoDenominator_==denomreco && numberOfValidTrackClusters != 0 ) quality=purity;
else quality=0;

if( quality>qualitySimToReco_ && !( threeHitTracksAreSpecial_ && numberOfSimulatedHits==3 && numberOfSharedHits<3.0 ) && ( absoluteNumberOfHits_ || (purity>puritySimToReco_) ) )
if( quality>qualitySimToReco_ && !( threeHitTracksAreSpecial_ && numberOfSimulatedHits==3 && numberOfSharedClusters<3.0 ) && ( absoluteNumberOfHits_ || (purity>puritySimToReco_) ) )
{
// Getting the RefToBase is dependent on the type of trackCollection, so delegate that to an overload.
returnValue.insert( trackingParticleRef, std::make_pair( ::getRefToTrackAt(trackCollection,i), quality ) );
Expand Down Expand Up @@ -516,6 +516,12 @@ template<typename iter> double QuickTrackAssociatorByHitsImpl::getDoubleCount( c
// This method is largely copied from the standard TrackAssociatorByHits. Once I've tested how much difference
// it makes I'll go through and comment it properly.

// FIXME: It may be that this piece is not fully correct for
// counting how many times a single *cluster* is matched to many
// SimTracks of a single TrackingParticle (see comments in
// getDoubleCount(ClusterTPAssociation) overload). To be verified
// some time.

double doubleCount=0.0;
std::vector < SimHitIdpr > SimTrackIdsDC;

Expand Down Expand Up @@ -552,16 +558,29 @@ template<typename iter> double QuickTrackAssociatorByHitsImpl::getDoubleCount( c
// This code here was written by Subir Sarkar. I'm just splitting it off into a
// separate method. - Grimes 01/May/2014

// The point here is that the electron TrackingParticles may contain
// multiple SimTracks (from the bremsstrahling), and (historically)
// the each matched hit/cluster has been multiplied by "how many
// SimTracks from the TrackingParticle" it contains charge from.
// Here the amount of this double counting is calculated, so that it
// can be subtracted by the calling code.
//
// Note that recently (hence "historically" in the paragraph above)
// the ClusterTPAssociationProducer was changed to remove the
// duplicate cluster->TP associations (hence making this function
// obsolete), but there is more recent proof that there is some
// duplication left (to be investigated).

double doubleCount=0;
std::vector < SimHitIdpr > SimTrackIdsDC;

for( iter iHit=startIterator; iHit != endIterator; iHit++ )
{
int idcount=0;

std::vector < OmniClusterRef > oClusters=getMatchedClusters( iHit, iHit + 1 ); //only for the cluster being checked
for( std::vector<OmniClusterRef>::const_iterator it=oClusters.begin(); it != oClusters.end(); ++it )
{
int idcount=0;

auto range = clusterToTPList.equal_range(*it);
if( range.first != range.second )
{
Expand All @@ -574,13 +593,13 @@ template<typename iter> double QuickTrackAssociatorByHitsImpl::getDoubleCount( c
}
}
}
}

if( idcount > 1 ) {
const auto subdetId = getHitFromIter(iHit)->geographicalId().subdetId();
const double weight = (subdetId == PixelSubdetector::PixelBarrel || subdetId == PixelSubdetector::PixelEndcap) ? pixelHitWeight_ : 1.0;
doubleCount += weight*(idcount - 1);
}
if( idcount > 1 ) {
const auto subdetId = getHitFromIter(iHit)->geographicalId().subdetId();
const double weight = (subdetId == PixelSubdetector::PixelBarrel || subdetId == PixelSubdetector::PixelEndcap) ? pixelHitWeight_ : 1.0;
doubleCount += weight*(idcount - 1);
}
}
}

return doubleCount;
Expand Down Expand Up @@ -609,24 +628,24 @@ reco::RecoToSimCollectionSeed QuickTrackAssociatorByHitsImpl::associateRecoToSim
++iTrackingParticleQualityPair )
{
const edm::Ref<TrackingParticleCollection>& trackingParticleRef=iTrackingParticleQualityPair->first;
double numberOfSharedHits=iTrackingParticleQualityPair->second;
double numberOfValidTrackHits=weightedNumberOfTrackHits(*pSeed);
double numberOfSharedClusters=iTrackingParticleQualityPair->second;
double numberOfValidTrackClusters = clusterToTPMap_ ? weightedNumberOfTrackClusters(*pSeed, *clusterToTPMap_) : weightedNumberOfTrackClusters(*pSeed, *hitAssociator_);

if( numberOfSharedHits == 0.0 ) continue; // No point in continuing if there was no association
if( numberOfSharedClusters == 0.0 ) continue; // No point in continuing if there was no association

//if electron subtract double counting
if( abs( trackingParticleRef->pdgId() ) == 11 && (trackingParticleRef->g4Track_end() - trackingParticleRef->g4Track_begin()) > 1 )
{
if( clusterToTPMap_ ) numberOfSharedHits-=getDoubleCount( *clusterToTPMap_, pSeed->recHits().first, pSeed->recHits().second, trackingParticleRef );
else numberOfSharedHits-=getDoubleCount( *hitAssociator_, pSeed->recHits().first, pSeed->recHits().second, trackingParticleRef );
if( clusterToTPMap_ ) numberOfSharedClusters-=getDoubleCount( *clusterToTPMap_, pSeed->recHits().first, pSeed->recHits().second, trackingParticleRef );
else numberOfSharedClusters-=getDoubleCount( *hitAssociator_, pSeed->recHits().first, pSeed->recHits().second, trackingParticleRef );
}

double quality;
if( absoluteNumberOfHits_ ) quality = numberOfSharedHits;
else if( numberOfValidTrackHits != 0.0 ) quality = numberOfSharedHits / numberOfValidTrackHits;
if( absoluteNumberOfHits_ ) quality = numberOfSharedClusters;
else if( numberOfValidTrackClusters != 0.0 ) quality = numberOfSharedClusters / numberOfValidTrackClusters;
else quality=0;

if( quality > cutRecoToSim_ && !(threeHitTracksAreSpecial_ && pSeed->nHits() == 3 && numberOfSharedHits < 3.0) )
if( quality > cutRecoToSim_ && !(threeHitTracksAreSpecial_ && pSeed->nHits() == 3 && numberOfSharedClusters < 3.0) )
{
returnValue.insert( edm::RefToBase < TrajectorySeed > (pSeedCollectionHandle_, i), std::make_pair( trackingParticleRef, quality ) );
}
Expand Down Expand Up @@ -668,20 +687,20 @@ reco::SimToRecoCollectionSeed QuickTrackAssociatorByHitsImpl::associateSimToReco
++iTrackingParticleQualityPair )
{
const edm::Ref<TrackingParticleCollection>& trackingParticleRef=iTrackingParticleQualityPair->first;
double numberOfSharedHits=iTrackingParticleQualityPair->second;
double numberOfValidTrackHits=weightedNumberOfTrackHits(*pSeed);
double numberOfSharedClusters=iTrackingParticleQualityPair->second;
double numberOfValidTrackClusters = clusterToTPMap_ ? weightedNumberOfTrackClusters(*pSeed, *clusterToTPMap_) :weightedNumberOfTrackClusters(*pSeed, *hitAssociator_);
size_t numberOfSimulatedHits=0; // Set a few lines below, but only if required.

if( numberOfSharedHits == 0.0 ) continue; // No point in continuing if there was no association
if( numberOfSharedClusters == 0.0 ) continue; // No point in continuing if there was no association

//if electron subtract double counting
if( abs( trackingParticleRef->pdgId() ) == 11 && (trackingParticleRef->g4Track_end() - trackingParticleRef->g4Track_begin()) > 1 )
{
if( clusterToTPMap_ ) numberOfSharedHits-=getDoubleCount( *clusterToTPMap_, pSeed->recHits().first, pSeed->recHits().second, trackingParticleRef );
else numberOfSharedHits-=getDoubleCount( *hitAssociator_, pSeed->recHits().first, pSeed->recHits().second, trackingParticleRef );
if( clusterToTPMap_ ) numberOfSharedClusters-=getDoubleCount( *clusterToTPMap_, pSeed->recHits().first, pSeed->recHits().second, trackingParticleRef );
else numberOfSharedClusters-=getDoubleCount( *hitAssociator_, pSeed->recHits().first, pSeed->recHits().second, trackingParticleRef );
}

if( simToRecoDenominator_ == denomsim || (numberOfSharedHits < 3.0 && threeHitTracksAreSpecial_) ) // the numberOfSimulatedHits is not always required, so can skip counting in some circumstances
if( simToRecoDenominator_ == denomsim || (numberOfSharedClusters < 3.0 && threeHitTracksAreSpecial_) ) // the numberOfSimulatedHits is not always required, so can skip counting in some circumstances
{
// Note that in the standard TrackAssociatorByHits, all of the hits in associatedTrackingParticleHits are checked for
// various things. I'm not sure what these checks are for but they depend on the UseGrouping and UseSplitting settings.
Expand All @@ -690,15 +709,15 @@ reco::SimToRecoCollectionSeed QuickTrackAssociatorByHitsImpl::associateSimToReco
numberOfSimulatedHits=trackingParticleRef->numberOfTrackerHits();
}

double purity = numberOfSharedHits / numberOfValidTrackHits;
double purity = numberOfSharedClusters / numberOfValidTrackClusters;
double quality;
if( absoluteNumberOfHits_ ) quality = numberOfSharedHits;
else if( simToRecoDenominator_ == denomsim && numberOfSimulatedHits != 0 ) quality= numberOfSharedHits
if( absoluteNumberOfHits_ ) quality = numberOfSharedClusters;
else if( simToRecoDenominator_ == denomsim && numberOfSimulatedHits != 0 ) quality= numberOfSharedClusters
/ static_cast<double>( numberOfSimulatedHits );
else if( simToRecoDenominator_ == denomreco && numberOfValidTrackHits != 0.0 ) quality=purity;
else if( simToRecoDenominator_ == denomreco && numberOfValidTrackClusters != 0.0 ) quality=purity;
else quality=0;

if( quality > qualitySimToReco_ && !(threeHitTracksAreSpecial_ && numberOfSimulatedHits == 3 && numberOfSharedHits < 3.0)
if( quality > qualitySimToReco_ && !(threeHitTracksAreSpecial_ && numberOfSimulatedHits == 3 && numberOfSharedClusters < 3.0)
&& (absoluteNumberOfHits_ || (purity > puritySimToReco_)) )
{
returnValue.insert( trackingParticleRef, std::make_pair( edm::RefToBase < TrajectorySeed > (pSeedCollectionHandle_, i), quality ) );
Expand All @@ -711,14 +730,15 @@ reco::SimToRecoCollectionSeed QuickTrackAssociatorByHitsImpl::associateSimToReco
return returnValue;
}

double QuickTrackAssociatorByHitsImpl::weightedNumberOfTrackHits(const reco::Track& track) const {
// count hits
double QuickTrackAssociatorByHitsImpl::weightedNumberOfTrackClusters(const reco::Track& track, const TrackerHitAssociator&) const {
const reco::HitPattern& p = track.hitPattern();
const auto pixelHits = p.numberOfValidPixelHits();
const auto otherHits = p.numberOfValidHits() - pixelHits;
return pixelHits*pixelHitWeight_ + otherHits;
}

double QuickTrackAssociatorByHitsImpl::weightedNumberOfTrackHits(const TrajectorySeed& seed) const {
double QuickTrackAssociatorByHitsImpl::weightedNumberOfTrackClusters(const TrajectorySeed& seed, const TrackerHitAssociator&) const {
double sum = 0.0;
for(auto iHit=seed.recHits().first; iHit!=seed.recHits().second; ++iHit) {
const auto subdetId = getHitFromIter(iHit)->geographicalId().subdetId();
Expand All @@ -727,3 +747,32 @@ double QuickTrackAssociatorByHitsImpl::weightedNumberOfTrackHits(const Trajector
}
return sum;
}

// count clusters
double QuickTrackAssociatorByHitsImpl::weightedNumberOfTrackClusters(const reco::Track& track, const ClusterTPAssociation&) const {
return weightedNumberOfTrackClusters(track.recHitsBegin(), track.recHitsEnd());
}
double QuickTrackAssociatorByHitsImpl::weightedNumberOfTrackClusters(const TrajectorySeed& seed, const ClusterTPAssociation&) const {
const auto& hitRange = seed.recHits();
return weightedNumberOfTrackClusters(hitRange.first, hitRange.second);
}

template<typename iter> double QuickTrackAssociatorByHitsImpl::weightedNumberOfTrackClusters(iter begin, iter end) const {

double weightedClusters = 0.0;
for (iter iRecHit = begin; iRecHit != end; ++iRecHit) {

const auto subdetId = getHitFromIter(iRecHit)->geographicalId().subdetId();
const double weight = (subdetId == PixelSubdetector::PixelBarrel || subdetId == PixelSubdetector::PixelEndcap) ? pixelHitWeight_ : 1.0;
LogTrace("QuickTrackAssociatorByHitsImpl") << " detId: " << getHitFromIter(iRecHit)->geographicalId().rawId();
LogTrace("QuickTrackAssociatorByHitsImpl") << " weight: " << weight;
std::vector < OmniClusterRef > oClusters=getMatchedClusters( iRecHit, iRecHit + 1 ); //only for the cluster being checked
for( std::vector<OmniClusterRef>::const_iterator it=oClusters.begin(); it != oClusters.end(); ++it ) {
weightedClusters += weight;
}
}
LogTrace("QuickTrackAssociatorByHitsImpl") << " total weighted clusters: " << weightedClusters;

return weightedClusters;
}

16 changes: 8 additions & 8 deletions SimTracker/TrackAssociatorProducers/plugins/QuickTrackAssociatorByHitsImpl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,26 +79,20 @@ class QuickTrackAssociatorByHitsImpl : public reco::TrackToTrackingParticleAssoc
bool threeHitTracksAreSpecial,
SimToRecoDenomType simToRecoDenominator);

virtual
reco::RecoToSimCollection associateRecoToSim( const edm::Handle<edm::View<reco::Track> >& trackCollectionHandle,
const edm::Handle<TrackingParticleCollection>& trackingParticleCollectionHandle) const override;
virtual
reco::SimToRecoCollection associateSimToReco( const edm::Handle<edm::View<reco::Track> >& trackCollectionHandle,
const edm::Handle<TrackingParticleCollection>& trackingParticleCollectionHandle) const override;
virtual
reco::RecoToSimCollection associateRecoToSim( const edm::RefToBaseVector<reco::Track>& trackCollection,
const edm::RefVector<TrackingParticleCollection>& trackingParticleCollection) const override;

virtual
reco::SimToRecoCollection associateSimToReco( const edm::RefToBaseVector<reco::Track>& trackCollection,
const edm::RefVector<TrackingParticleCollection>& trackingParticleCollection) const override;

//seed
virtual
reco::RecoToSimCollectionSeed associateRecoToSim(const edm::Handle<edm::View<TrajectorySeed> >&,
const edm::Handle<TrackingParticleCollection>&) const override;

virtual
reco::SimToRecoCollectionSeed associateSimToReco(const edm::Handle<edm::View<TrajectorySeed> >&,
const edm::Handle<TrackingParticleCollection>&) const override;

Expand Down Expand Up @@ -177,8 +171,14 @@ class QuickTrackAssociatorByHitsImpl : public reco::TrackToTrackingParticleAssoc
return &(*iter);
}

double weightedNumberOfTrackHits(const reco::Track& track) const;
double weightedNumberOfTrackHits(const TrajectorySeed& seed) const;
// The last parameter is used to decide whether we cound hits or clusters
double weightedNumberOfTrackClusters(const reco::Track& track, const TrackerHitAssociator&) const;
double weightedNumberOfTrackClusters(const TrajectorySeed& seed, const TrackerHitAssociator&) const;
double weightedNumberOfTrackClusters(const reco::Track& track, const ClusterTPAssociation&) const;
double weightedNumberOfTrackClusters(const TrajectorySeed& seed, const ClusterTPAssociation&) const;

// called only by weightedNumberOfTrackClusters(..., ClusterTPAssociation)
template<typename iter> double weightedNumberOfTrackClusters(iter begin, iter end) const ;

/** @brief creates either a ClusterTPAssociation OR a TrackerHitAssociator and stores it in the provided unique_ptr. The other will be null.
*
Expand Down