FastTSGFromIOHit.cc

Go to the documentation of this file.

#include "FastSimulation/Muons/plugins/FastTSGFromIOHit.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Framework/interface/Event.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h"
#include "DataFormats/MuonSeed/interface/L3MuonTrajectorySeed.h"
#include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h"
#include "DataFormats/Math/interface/deltaPhi.h"
#include "FastSimulation/Tracking/interface/FastTrackingHelper.h"

FastTSGFromIOHit::FastTSGFromIOHit(const edm::ParameterSet & iConfig,edm::ConsumesCollector& iC) 
{

  theCategory = "FastSimulation|Muons||FastTSGFromIOHit";

  thePtCut = iConfig.getParameter<double>("PtCut");

  simTracksTk = iC.consumes<edm::SimTrackContainer>(iConfig.getParameter<edm::InputTag>("SimTrackCollectionLabel"));
  const auto & seedLabels = iConfig.getParameter<std::vector<edm::InputTag> >("SeedCollectionLabels");
  for(const auto & seedLabel : seedLabels)
  {
      seedsTks.push_back(iC.consumes<TrajectorySeedCollection>(seedLabel));
  }

}

FastTSGFromIOHit::~FastTSGFromIOHit()
{

  LogTrace(theCategory) << " FastTSGFromIOHit dtor called ";

}

void FastTSGFromIOHit::trackerSeeds(const TrackCand& staMuon, const TrackingRegion& region, const TrackerTopology *tTopo, std::vector<TrajectorySeed> & result) 
{
    // Make a ref to l2 muon
    reco::TrackRef muRef(staMuon.second);

    // Cut on muons with low momenta
    if ( muRef->pt() < thePtCut 
     || muRef->innerMomentum().Rho() < thePtCut 
     || muRef->innerMomentum().R() < 2.5 )
    {
    return;
    }
    
    // Retrieve the Monte Carlo truth (SimTracks)
    edm::Handle<edm::SimTrackContainer> simTracks;
    getEvent()->getByToken(simTracksTk,simTracks);
    
    // Retrieve Seed collection
    std::vector<edm::Handle<TrajectorySeedCollection> > seedCollections;
    seedCollections.resize(seedsTks.size());
    for ( unsigned iSeed = 0 ; iSeed < seedsTks.size() ; iSeed++) 
    {
    getEvent()->getByToken(seedsTks[iSeed],seedCollections[iSeed]);
    }
    
    // cast the tracking region
    const RectangularEtaPhiTrackingRegion & regionRef = dynamic_cast<const RectangularEtaPhiTrackingRegion & > (region);
  
    // select and store seeds
    std::set<unsigned> simTrackIds;
    for ( const auto & seeds : seedCollections) 
    { 
    for ( const auto & seed : *seeds)
    {
        // Find the simtrack corresponding to the seed
        TrajectorySeed::range recHitRange = seed.recHits();
        const FastTrackerRecHit * firstRecHit = (const FastTrackerRecHit*) (&(*(recHitRange.first)));
        int simTrackId = firstRecHit->simTrackId(0);
        const SimTrack & simTrack = (*simTracks)[simTrackId];
      
        // skip if simTrack already associated to a seed
        if( simTrackIds.find(simTrackId) != simTrackIds.end() )
        {
        continue;   
        }
        simTrackIds.insert(simTrackId);
      
        // filter seed
        if( !clean(muRef,regionRef,&seed,simTrack) )
        {
        continue;
        }

        // store results
        result.push_back(L3MuonTrajectorySeed(seed,muRef));
    } 
    } 
}

bool
FastTSGFromIOHit::clean(reco::TrackRef muRef,
            const RectangularEtaPhiTrackingRegion& region,
            const TrajectorySeed* aSeed,
            const SimTrack& theSimTrack) 
{

  // Eta cleaner   
  const PixelRecoRange<float>& etaRange = region.etaRange() ;  

  double etaSeed = theSimTrack.momentum().Eta();
  double etaLimit  = (fabs(fabs(etaRange.max())-fabs(etaRange.mean())) <0.05) ? 
    0.05 : fabs(fabs(etaRange.max()) - fabs(etaRange.mean())) ;
  bool inEtaRange = 
    etaSeed >= (etaRange.mean() - etaLimit) && 
    etaSeed <= (etaRange.mean() + etaLimit) ;
  if  ( !inEtaRange )
  {
      return false;
  }
  
  // Phi cleaner
  const TkTrackingRegionsMargin<float>& phiMargin = region.phiMargin();
  double phiSeed = theSimTrack.momentum().Phi(); 
  double phiLimit  = (phiMargin.right() < 0.05 ) ? 0.05 : phiMargin.right(); 
  bool inPhiRange = 
    (fabs(deltaPhi(phiSeed,double(region.direction().phi()))) < phiLimit );
  if  ( !inPhiRange )
  {
      return false;
  }
  
  // pt cleaner
  double ptSeed  = std::sqrt(theSimTrack.momentum().Perp2());
  double ptMin   = (region.ptMin()>3.5) ? 3.5: region.ptMin();  
  bool inPtRange = ptSeed >= ptMin &&  ptSeed<= 2*(muRef->pt());
  if  ( !inPtRange )
  {
      return false;
  }
  return true;
  
}