/data/refman/pasoursint/CMSSW_5_3_3/src/SimTracker/TrackAssociation/src/QuickTrackAssociatorByHits.cc

Go to the documentation of this file.

#include "SimTracker/TrackAssociation/interface/QuickTrackAssociatorByHits.h"

#include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

QuickTrackAssociatorByHits::QuickTrackAssociatorByHits( const edm::ParameterSet& config )
        : pHitAssociator_(NULL), pEventForWhichAssociatorIsValid_(NULL),
          absoluteNumberOfHits_( config.getParameter<bool>( "AbsoluteNumberOfHits" ) ),
          qualitySimToReco_( config.getParameter<double>( "Quality_SimToReco" ) ),
          puritySimToReco_( config.getParameter<double>( "Purity_SimToReco" ) ),
          cutRecoToSim_( config.getParameter<double>( "Cut_RecoToSim" ) ),
          threeHitTracksAreSpecial_( config.getParameter<bool> ( "ThreeHitTracksAreSpecial" ) )
{
        //
        // Check whether the denominator when working out the percentage of shared hits should
        // be the number of simulated hits or the number of reconstructed hits.
        //
        std::string denominatorString=config.getParameter<std::string>("SimToRecoDenominator");
        if( denominatorString=="sim" ) simToRecoDenominator_=denomsim;
        else if( denominatorString=="reco" ) simToRecoDenominator_=denomreco;
        else throw cms::Exception( "QuickTrackAssociatorByHits" ) << "SimToRecoDenominator not specified as sim or reco";

        //
        // Set up the parameter set for the hit associator
        //
        hitAssociatorParameters_.addParameter<bool>( "associatePixel", config.getParameter<bool>("associatePixel") );
        hitAssociatorParameters_.addParameter<bool>( "associateStrip", config.getParameter<bool>("associateStrip") );
        // This is the important one, it stops the hit associator searching through the list of sim hits.
        // I only want to use the hit associator methods that work on the hit IDs (i.e. the uint32_t trackId
        // and the EncodedEventId eventId) so I'm not interested in matching that to the PSimHit objects.
        hitAssociatorParameters_.addParameter<bool>("associateRecoTracks",true);

        //
        // Do some checks on whether UseGrouped or UseSplitting have been set. They're not used
        // unlike the standard TrackAssociatorByHits so show a warning.
        //
        bool useGrouped, useSplitting;
        if( config.exists("UseGrouped") ) useGrouped=config.getParameter<bool>("UseGrouped");
        else useGrouped=true;

        if( config.exists("UseSplitting") ) useSplitting=config.getParameter<bool>("UseSplitting");
        else useSplitting=true;

        // This associator works as though both UseGrouped and UseSplitting were set to true, so show a
        // warning if this isn't the case.
        if( !(useGrouped && useSplitting) )
        {
                edm::LogWarning("QuickTrackAssociatorByHits") << "UseGrouped and/or UseSplitting has been set to false, but this associator ignores that setting.";
        }
}

QuickTrackAssociatorByHits::~QuickTrackAssociatorByHits()
{
        delete pHitAssociator_;
}

QuickTrackAssociatorByHits::QuickTrackAssociatorByHits( const QuickTrackAssociatorByHits& otherAssociator )
        : pEventForWhichAssociatorIsValid_(otherAssociator.pEventForWhichAssociatorIsValid_),
          hitAssociatorParameters_(otherAssociator.hitAssociatorParameters_),
          absoluteNumberOfHits_(otherAssociator.absoluteNumberOfHits_),
          qualitySimToReco_(otherAssociator.qualitySimToReco_),
          puritySimToReco_(otherAssociator.puritySimToReco_),
          cutRecoToSim_(otherAssociator.cutRecoToSim_),
          threeHitTracksAreSpecial_(otherAssociator.threeHitTracksAreSpecial_),
          simToRecoDenominator_(otherAssociator.simToRecoDenominator_),
          pTrackCollectionHandle_(otherAssociator.pTrackCollectionHandle_),
          pTrackCollection_(otherAssociator.pTrackCollection_),
          pTrackingParticleCollectionHandle_(otherAssociator.pTrackingParticleCollectionHandle_),
          pTrackingParticleCollection_(otherAssociator.pTrackingParticleCollection_)

{
        // No operation other than the initialiser list. That copies everything straight from the other
        // associator, except for pHitAssociator_ which needs a deep copy or both instances will try
        // and free it on deletion.  If it wasn't for pHitAssociator_ the default copy constructor and
        // assignment operator would be sufficient.

        // Actually, need to check the other hit associator isn't null or the pointer dereference would
        // probably cause a segmentation fault.
        if( otherAssociator.pHitAssociator_ ) pHitAssociator_=new TrackerHitAssociator(*otherAssociator.pHitAssociator_);
        else pHitAssociator_=NULL;
}

QuickTrackAssociatorByHits& QuickTrackAssociatorByHits::operator=( const QuickTrackAssociatorByHits& otherAssociator )
{
        // Free up the old pHitAssociator_
        delete pHitAssociator_;

        //
        // pHitAssociator_ needs to be given a deep copy of the object, but everything else can
        // can be shallow copied from the other associator.
        //
        if( otherAssociator.pHitAssociator_ ) pHitAssociator_=new TrackerHitAssociator(*otherAssociator.pHitAssociator_);
        else pHitAssociator_=NULL;
        pEventForWhichAssociatorIsValid_=otherAssociator.pEventForWhichAssociatorIsValid_;
        hitAssociatorParameters_=otherAssociator.hitAssociatorParameters_;
        absoluteNumberOfHits_=otherAssociator.absoluteNumberOfHits_;
        qualitySimToReco_=otherAssociator.qualitySimToReco_;
        puritySimToReco_=otherAssociator.puritySimToReco_;
        cutRecoToSim_=otherAssociator.cutRecoToSim_;
        threeHitTracksAreSpecial_=otherAssociator.threeHitTracksAreSpecial_;
        simToRecoDenominator_=otherAssociator.simToRecoDenominator_;
        pTrackCollectionHandle_=otherAssociator.pTrackCollectionHandle_;
        pTrackCollection_=otherAssociator.pTrackCollection_;
        pTrackingParticleCollectionHandle_=otherAssociator.pTrackingParticleCollectionHandle_;
        pTrackingParticleCollection_=otherAssociator.pTrackingParticleCollection_;

        return *this;
}

reco::RecoToSimCollection QuickTrackAssociatorByHits::associateRecoToSim( edm::Handle<edm::View<reco::Track> >& trackCollectionHandle,
                                                                              edm::Handle<TrackingParticleCollection>& trackingParticleCollectionHandle,
                                                                              const edm::Event* pEvent,
                                                                              const edm::EventSetup* pSetup ) const
{
        initialiseHitAssociator( pEvent );
        pTrackCollectionHandle_=&trackCollectionHandle;
        pTrackingParticleCollectionHandle_=&trackingParticleCollectionHandle;
        pTrackCollection_=NULL;
        pTrackingParticleCollection_=NULL;

        // This method checks which collection type is set to NULL, and uses the other one.
        return associateRecoToSimImplementation();
}

reco::SimToRecoCollection QuickTrackAssociatorByHits::associateSimToReco( edm::Handle<edm::View<reco::Track> >& trackCollectionHandle,
                                                                              edm::Handle<TrackingParticleCollection>& trackingParticleCollectionHandle,
                                                                              const edm::Event * pEvent,
                                                                              const edm::EventSetup * pSetup ) const
{
        initialiseHitAssociator( pEvent );
        pTrackCollectionHandle_=&trackCollectionHandle;
        pTrackingParticleCollectionHandle_=&trackingParticleCollectionHandle;
        pTrackCollection_=NULL;
        pTrackingParticleCollection_=NULL;

        // This method checks which collection type is set to NULL, and uses the other one.
        return associateSimToRecoImplementation();
}


reco::RecoToSimCollection QuickTrackAssociatorByHits::associateRecoToSim(const edm::RefToBaseVector<reco::Track>& trackCollection,
                                                                             const edm::RefVector<TrackingParticleCollection>& trackingParticleCollection,
                                                                             const edm::Event* pEvent,
                                                                             const edm::EventSetup* pSetup ) const
{
        initialiseHitAssociator( pEvent );
        pTrackCollectionHandle_=NULL;
        pTrackingParticleCollectionHandle_=NULL;
        pTrackCollection_=&trackCollection;
        pTrackingParticleCollection_=&trackingParticleCollection;

        // This method checks which collection type is set to NULL, and uses the other one.
        return associateRecoToSimImplementation();
}

reco::SimToRecoCollection QuickTrackAssociatorByHits::associateSimToReco(const edm::RefToBaseVector<reco::Track>& trackCollection,
                                                                             const edm::RefVector<TrackingParticleCollection>& trackingParticleCollection,
                                                                             const edm::Event* pEvent,
                                                                             const edm::EventSetup* pSetup ) const
{
        initialiseHitAssociator( pEvent );
        pTrackCollectionHandle_=NULL;
        pTrackingParticleCollectionHandle_=NULL;
        pTrackCollection_=&trackCollection;
        pTrackingParticleCollection_=&trackingParticleCollection;

        // This method checks which collection type is set to NULL, and uses the other one.
        return associateSimToRecoImplementation();
}

reco::RecoToSimCollection QuickTrackAssociatorByHits::associateRecoToSimImplementation() const
{
        reco::RecoToSimCollection returnValue;

        size_t collectionSize;
        // Need to check which pointer is valid to get the collection size
        if( pTrackCollection_ ) collectionSize=pTrackCollection_->size();
        else collectionSize=(*pTrackCollectionHandle_)->size();

        for( size_t i=0; i<collectionSize; ++i )
        {
                const reco::Track* pTrack; // Get a normal pointer for ease of use.
                if( pTrackCollection_ ) pTrack=&*(*pTrackCollection_)[i]; // Possibly the most obscure dereference I've ever had to write
                else pTrack=&(*pTrackCollectionHandle_->product())[i];

                // The return of this function has first as the index and second as the number of associated hits
                std::vector< std::pair<edm::Ref<TrackingParticleCollection>,size_t> > trackingParticleQualityPairs=associateTrack( pTrack );
                for( std::vector< std::pair<edm::Ref<TrackingParticleCollection>,size_t> >::const_iterator iTrackingParticleQualityPair=trackingParticleQualityPairs.begin();
                                                iTrackingParticleQualityPair!=trackingParticleQualityPairs.end(); ++iTrackingParticleQualityPair )
                {
                        const edm::Ref<TrackingParticleCollection>& trackingParticleRef=iTrackingParticleQualityPair->first;
                        size_t numberOfSharedHits=iTrackingParticleQualityPair->second;
                        size_t numberOfValidTrackHits=pTrack->found();

                        if( numberOfSharedHits==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( pTrack->recHitsBegin(), pTrack->recHitsEnd(), *trackingParticleRef );
                        }

                        double quality;
                        if( absoluteNumberOfHits_ ) quality=static_cast<double>( numberOfSharedHits );
                        else if( numberOfValidTrackHits != 0 ) quality=(static_cast<double>(numberOfSharedHits) / static_cast<double>(numberOfValidTrackHits) );
                        else quality=0;

                        if( quality > cutRecoToSim_ && !( threeHitTracksAreSpecial_ && numberOfValidTrackHits==3 && numberOfSharedHits<3 ) )
                        {
                                if( pTrackCollection_ ) returnValue.insert( (*pTrackCollection_)[i], std::make_pair( trackingParticleRef, quality ));
                                else returnValue.insert( edm::RefToBase<reco::Track>(*pTrackCollectionHandle_,i), std::make_pair( trackingParticleRef, quality ));
                        }
                }
        }
        return returnValue;
}

reco::SimToRecoCollection QuickTrackAssociatorByHits::associateSimToRecoImplementation() const
{
        reco::SimToRecoCollection returnValue;

        size_t collectionSize;
        // Need to check which pointer is valid to get the collection size
        if( pTrackCollection_ ) collectionSize=pTrackCollection_->size();
        else collectionSize=(*pTrackCollectionHandle_)->size();

        for( size_t i=0; i<collectionSize; ++i )
        {
                const reco::Track* pTrack; // Get a normal pointer for ease of use.
                if( pTrackCollection_ ) pTrack=&*(*pTrackCollection_)[i]; // Possibly the most obscure dereference I've ever had to write
                else pTrack=&(*pTrackCollectionHandle_->product())[i];

                // The return of this function has first as an edm:Ref to the associated TrackingParticle, and second as the number of associated hits
                std::vector< std::pair<edm::Ref<TrackingParticleCollection>,size_t> > trackingParticleQualityPairs=associateTrack( pTrack );
                for( std::vector< std::pair<edm::Ref<TrackingParticleCollection>,size_t> >::const_iterator iTrackingParticleQualityPair=trackingParticleQualityPairs.begin();
                                iTrackingParticleQualityPair!=trackingParticleQualityPairs.end(); ++iTrackingParticleQualityPair )
                {
                        const edm::Ref<TrackingParticleCollection>& trackingParticleRef=iTrackingParticleQualityPair->first;
                        size_t numberOfSharedHits=iTrackingParticleQualityPair->second;
                        size_t numberOfValidTrackHits=pTrack->found();
                        size_t numberOfSimulatedHits=0; // Set a few lines below, but only if required.

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

                        if( simToRecoDenominator_==denomsim || (numberOfSharedHits<3 && 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.
                                // This associator works as though both UseGrouping and UseSplitting were set to true, i.e. just counts the number of
                                // hits in the tracker.
                                numberOfSimulatedHits=trackingParticleRef->trackPSimHit(DetId::Tracker).size();
                        }

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

                        if( quality>qualitySimToReco_ && !( threeHitTracksAreSpecial_ && numberOfSimulatedHits==3 && numberOfSharedHits<3 ) && ( absoluteNumberOfHits_ || (purity>puritySimToReco_) ) )
                        {
                                if( pTrackCollection_ ) returnValue.insert( trackingParticleRef, std::make_pair( (*pTrackCollection_)[i], quality ) );
                                else returnValue.insert( trackingParticleRef, std::make_pair( edm::RefToBase<reco::Track>(*pTrackCollectionHandle_,i) , quality ) );
                        }
                }
        }
        return returnValue;

}

std::vector< std::pair<edm::Ref<TrackingParticleCollection>,size_t> > QuickTrackAssociatorByHits::associateTrack( const reco::Track* pTrack ) const
{
        // The pairs in this vector have a Ref to the associated TrackingParticle as "first" and the number of associated hits as "second"
        std::vector< std::pair<edm::Ref<TrackingParticleCollection>,size_t> > returnValue;

        // The pairs in this vector have first as the sim track identifiers, and second the number of reco hits associated to that sim track.
        // 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
        // number of reco hits.  The pair::second entries should add up to the total number of reco hits though.
        std::vector< std::pair<SimTrackIdentifiers,size_t> > hitIdentifiers=getAllSimTrackIdentifiers(pTrack);

        // Loop over the TrackingParticles
        size_t collectionSize;
        if( pTrackingParticleCollection_ ) collectionSize=pTrackingParticleCollection_->size();
        else collectionSize=(*pTrackingParticleCollectionHandle_)->size();

        for( size_t i=0; i<collectionSize; ++i )
        {
                const TrackingParticle* pTrackingParticle; // Convert to raw pointer for ease of use
                if( pTrackingParticleCollection_ ) pTrackingParticle=&*(*pTrackingParticleCollection_)[i];
                else pTrackingParticle=&(*pTrackingParticleCollectionHandle_->product())[i];

                // Ignore TrackingParticles with no hits
                if( pTrackingParticle->trackPSimHit().empty() ) continue;

                size_t numberOfAssociatedHits=0;
                // Loop over all of the sim track identifiers and see if any of them are part of this TrackingParticle. If they are, add
                // the number of reco hits associated to that sim track to the total number of associated hits.
                for( std::vector< std::pair<SimTrackIdentifiers,size_t> >::const_iterator iIdentifierCountPair=hitIdentifiers.begin(); iIdentifierCountPair!=hitIdentifiers.end(); ++iIdentifierCountPair )
                {
                        if( trackingParticleContainsIdentifier( pTrackingParticle, iIdentifierCountPair->first ) ) numberOfAssociatedHits+=iIdentifierCountPair->second;
                }

                if( numberOfAssociatedHits>0 )
                {
                        if( pTrackingParticleCollection_ ) returnValue.push_back( std::make_pair( (*pTrackingParticleCollection_)[i], numberOfAssociatedHits ) );
                        else returnValue.push_back( std::make_pair( edm::Ref<TrackingParticleCollection>( *pTrackingParticleCollectionHandle_, i ), numberOfAssociatedHits ) );
                }
        }

        return returnValue;
}

std::vector< std::pair<QuickTrackAssociatorByHits::SimTrackIdentifiers,size_t> > QuickTrackAssociatorByHits::getAllSimTrackIdentifiers( const reco::Track* pTrack ) const
{
        // The pairs in this vector have first as the sim track identifiers, and second the number of reco hits associated to that sim track.
        std::vector< std::pair<SimTrackIdentifiers,size_t> > returnValue;

        std::vector<SimTrackIdentifiers> simTrackIdentifiers;
        // Loop over all of the rec hits in the track
        for( trackingRecHit_iterator iRecHit=pTrack->recHitsBegin(); iRecHit!=pTrack->recHitsEnd(); ++iRecHit )
        {
                if( (*iRecHit)->isValid() )
                {
                        simTrackIdentifiers.clear();

                        // Get the identifiers for the sim track that this hit came from. There should only be one entry unless clusters
                        // have merged (as far as I know).
                        pHitAssociator_->associateHitId( **iRecHit, simTrackIdentifiers ); // This call fills simTrackIdentifiers

                        // Loop over each identifier, and add it to the return value only if it's not already in there
                        for( std::vector<SimTrackIdentifiers>::const_iterator iIdentifier=simTrackIdentifiers.begin(); iIdentifier!=simTrackIdentifiers.end(); ++iIdentifier )
                        {
                                std::vector< std::pair<SimTrackIdentifiers,size_t> >::iterator iIdentifierCountPair;
                                for( iIdentifierCountPair=returnValue.begin(); iIdentifierCountPair!=returnValue.end(); ++iIdentifierCountPair )
                                {
                                        if( iIdentifierCountPair->first.first==iIdentifier->first && iIdentifierCountPair->first.second==iIdentifier->second )
                                        {
                                                // This sim track identifier is already in the list, so increment the count of how many hits it relates to.
                                                ++iIdentifierCountPair->second;
                                                break;
                                        }
                                }
                                if( iIdentifierCountPair==returnValue.end() ) returnValue.push_back( std::make_pair(*iIdentifier,1) ); // This identifier wasn't found, so add it
                        }
                }
        }

        return returnValue;
}

bool QuickTrackAssociatorByHits::trackingParticleContainsIdentifier( const TrackingParticle* pTrackingParticle, const SimTrackIdentifiers& identifier ) const
{
        // Loop over all of the g4 tracks in the tracking particle
        for( std::vector<SimTrack>::const_iterator iSimTrack=pTrackingParticle->g4Track_begin(); iSimTrack!=pTrackingParticle->g4Track_end(); ++iSimTrack )
        {
                // And see if the sim track identifiers match
                if( iSimTrack->eventId()==identifier.second && iSimTrack->trackId()==identifier.first )
                {
                        return true;
                }
        }

        // If control has made it this far then none of the identifiers were found in
        // any of the g4 tracks, so return false.
        return false;
}

int QuickTrackAssociatorByHits::getDoubleCount( trackingRecHit_iterator startIterator, trackingRecHit_iterator endIterator, const TrackingParticle& associatedTrackingParticle ) const
{
        // 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.

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

        for( trackingRecHit_iterator iHit=startIterator; iHit != endIterator; iHit++ )
        {
                int idcount=0;
                SimTrackIdsDC.clear();
                pHitAssociator_->associateHitId( **iHit, SimTrackIdsDC );

                if( SimTrackIdsDC.size() > 1 )
                {
                        for( TrackingParticle::g4t_iterator g4T=associatedTrackingParticle.g4Track_begin(); g4T != associatedTrackingParticle.g4Track_end(); ++g4T )
                        {
                                if( find( SimTrackIdsDC.begin(), SimTrackIdsDC.end(), SimHitIdpr( ( *g4T).trackId(), SimTrackIdsDC.begin()->second ) )
                                                != SimTrackIdsDC.end() )
                                {
                                        idcount++;
                                }
                        }
                }
                if( idcount > 1 ) doubleCount+=(idcount - 1);
        }

        return doubleCount;
}


void QuickTrackAssociatorByHits::initialiseHitAssociator( const edm::Event* pEvent ) const
{
        // The intention of this function was to check whether the hit associator is still valid
        // (since in general associateSimToReco and associateRecoToSim are called for the same
        // event). I was doing this by recording the event pointer and checking it hasn't changed
        // but this doesn't appear to work. Until I find a way of uniquely identifying an event
        // I'll just create it anew each time.
//      if( pEventForWhichAssociatorIsValid_==pEvent && pEventForWhichAssociatorIsValid_!=NULL ) return; // Already set up so no need to do anything

        // Free up the previous instantiation
        delete pHitAssociator_;

        // Create a new instantiation using the new event
        pHitAssociator_=new TrackerHitAssociator( *pEvent, hitAssociatorParameters_ );
        pEventForWhichAssociatorIsValid_=pEvent;
}