GsfElectronTools.cc

Go to the documentation of this file.

#include "DataFormats/EgammaCandidates/interface/GsfElectronCoreFwd.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectronCore.h"
#include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrack.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
#include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
#include "DataFormats/EgammaReco/interface/ElectronSeed.h"
#include "RecoEgamma/EgammaElectronAlgos/interface/GsfElectronTools.h"

namespace gsfElectronTools {

  using namespace reco ;

  //=======================================================================================
  // Code from Puneeth Kalavase
  //=======================================================================================

  std::pair<TrackRef,float> getClosestCtfToGsf( GsfTrackRef const& gsfTrackRef,
                                                edm::Handle<reco::TrackCollection> const& ctfTracksH )
  {
    float maxFracShared = 0;
    TrackRef ctfTrackRef = TrackRef() ;
    const TrackCollection * ctfTrackCollection = ctfTracksH.product() ;

    // get the Hit Pattern for the gsfTrack
    const HitPattern& gsfHitPattern = gsfTrackRef->hitPattern();

    unsigned int counter = 0;
    for (auto  ctfTkIter = ctfTrackCollection->begin();
          ctfTkIter != ctfTrackCollection->end() ; ctfTkIter++, counter++ )
     {

      double dEta = gsfTrackRef->eta() - ctfTkIter->eta();
      double dPhi = gsfTrackRef->phi() - ctfTkIter->phi();
      double pi = acos(-1.);
      if(std::abs(dPhi) > pi) dPhi = 2*pi - std::abs(dPhi);

      // dont want to look at every single track in the event!
      if(sqrt(dEta*dEta + dPhi*dPhi) > 0.3) continue;

      unsigned int shared = 0 ;
      int gsfHitCounter = 0 ;
      int numGsfInnerHits = 0 ;
      int numCtfInnerHits = 0 ;
      // get the CTF Track Hit Pattern
      const HitPattern& ctfHitPattern = ctfTkIter->hitPattern();

      for ( auto elHitsIt = gsfTrackRef->recHitsBegin() ;
            elHitsIt != gsfTrackRef->recHitsEnd() ;
            elHitsIt++, gsfHitCounter++ )
       {
        if(!((**elHitsIt).isValid()))  //count only valid Hits
         { continue ; }

        // look only in the pixels/TIB/TID
        uint32_t gsfHit = gsfHitPattern.getHitPattern(HitPattern::TRACK_HITS, gsfHitCounter);
        if (!(HitPattern::pixelHitFilter(gsfHit) ||
            HitPattern::stripTIBHitFilter(gsfHit) ||
            HitPattern::stripTIDHitFilter(gsfHit)))
         { continue ; }

        numGsfInnerHits++ ;

        int ctfHitsCounter = 0 ;
        numCtfInnerHits = 0 ;
        for ( auto ctfHitsIt = ctfTkIter->recHitsBegin() ;
              ctfHitsIt != ctfTkIter->recHitsEnd() ;
              ctfHitsIt++, ctfHitsCounter++ )
         {
          if(!((**ctfHitsIt).isValid())) //count only valid Hits!
           { continue ; }

        uint32_t ctfHit = ctfHitPattern.getHitPattern(HitPattern::TRACK_HITS, ctfHitsCounter);
        if(!(HitPattern::pixelHitFilter(ctfHit) ||
              HitPattern::stripTIBHitFilter(ctfHit) ||
              HitPattern::stripTIDHitFilter(ctfHit)))
         { continue ; }

        numCtfInnerHits++ ;

          if( (**elHitsIt).sharesInput(&(**ctfHitsIt),TrackingRecHit::all) )
           {
            shared++ ;
            break ;
           }

         } //ctfHits iterator

       } //gsfHits iterator

      if ((numGsfInnerHits==0)||(numCtfInnerHits==0))
       { continue ; }

      if ( static_cast<float>(shared)/std::min(numGsfInnerHits,numCtfInnerHits) > maxFracShared )
       {
        maxFracShared = static_cast<float>(shared)/std::min(numGsfInnerHits, numCtfInnerHits);
        ctfTrackRef = TrackRef(ctfTracksH,counter);
       }

     } //ctfTrack iterator

    return make_pair(ctfTrackRef,maxFracShared) ;
  }

}