#include <FastTSGFromL2Muon.h>

Inheritance diagram for FastTSGFromL2Muon:

Public Member Functions
virtual void	beginRun (edm::Run &run, edm::EventSetup const &es)
	FastTSGFromL2Muon (const edm::ParameterSet &cfg)
virtual void	produce (edm::Event &ev, const edm::EventSetup &es)
virtual	~FastTSGFromL2Muon ()
Private Member Functions
bool	clean (reco::TrackRef muRef, RectangularEtaPhiTrackingRegion region, const BasicTrajectorySeed aSeed, const SimTrack &theSimTrack)
Private Attributes
edm::ParameterSet	theConfig
edm::InputTag	theL2CollectionLabel
double	thePtCut
MuonTrackingRegionBuilder *	theRegionBuilder
std::vector< edm::InputTag >	theSeedCollectionLabels
MuonServiceProxy *	theService
edm::InputTag	theSimTrackCollectionLabel

Detailed Description

Definition at line 22 of file FastTSGFromL2Muon.h.

Constructor & Destructor Documentation

FastTSGFromL2Muon::FastTSGFromL2Muon ( const edm::ParameterSet & cfg )

Definition at line 21 of file FastTSGFromL2Muon.cc.

References edm::ParameterSet::getParameter(), MuonServiceProxy_cff::MuonServiceProxy, theL2CollectionLabel, thePtCut, theSeedCollectionLabels, theService, and theSimTrackCollectionLabel.

                                                               : theConfig(cfg) 
{
  produces<L3MuonTrajectorySeedCollection>();

  edm::ParameterSet serviceParameters = 
    cfg.getParameter<edm::ParameterSet>("ServiceParameters");
  theService = new MuonServiceProxy(serviceParameters);

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

  theL2CollectionLabel = cfg.getParameter<edm::InputTag>("MuonCollectionLabel");
  theSeedCollectionLabels = cfg.getParameter<std::vector<edm::InputTag> >("SeedCollectionLabels");
  theSimTrackCollectionLabel  = cfg.getParameter<edm::InputTag>("SimTrackCollectionLabel");
  // useTFileService_ = cfg.getUntrackedParameter<bool>("UseTFileService",false);

}

FastTSGFromL2Muon::~FastTSGFromL2Muon ( ) [virtual]

Definition at line 38 of file FastTSGFromL2Muon.cc.

{
}

Member Function Documentation

void FastTSGFromL2Muon::beginRun	(	edm::Run &	run,
		edm::EventSetup const &	es
	)		`[virtual]`

Reimplemented from edm::EDProducer.

Definition at line 43 of file FastTSGFromL2Muon.cc.

References edm::ParameterSet::getParameter(), theConfig, theRegionBuilder, theService, and MuonServiceProxy::update().

{
  //update muon proxy service
  theService->update(es);
  
  //region builder
  edm::ParameterSet regionBuilderPSet = 
    theConfig.getParameter<edm::ParameterSet>("MuonTrackingRegionBuilder");
  theRegionBuilder = new MuonTrackingRegionBuilder(regionBuilderPSet,theService);
  
  /*
  if(useTFileService_) {
    edm::Service<TFileService> fs;
    h_nSeedPerTrack = fs->make<TH1F>("nSeedPerTrack","nSeedPerTrack",76,-0.5,75.5);
    h_nGoodSeedPerTrack = fs->make<TH1F>("nGoodSeedPerTrack","nGoodSeedPerTrack",75,-0.5,75.5);
    h_nGoodSeedPerEvent = fs->make<TH1F>("nGoodSeedPerEvent","nGoodSeedPerEvent",75,-0.5,75.5);
  } else {
    h_nSeedPerTrack = 0;
    h_nGoodSeedPerEvent = 0;
    h_nGoodSeedPerTrack = 0;
  }
  */

}

bool FastTSGFromL2Muon::clean	(	reco::TrackRef	muRef,
		RectangularEtaPhiTrackingRegion *	region,
		const BasicTrajectorySeed *	aSeed,
		const SimTrack &	theSimTrack
	)		`[private]`

Definition at line 180 of file FastTSGFromL2Muon.cc.

References Geom::deltaPhi(), TrackingRegionBase::direction(), MuonErrorMatrixValues_cff::etaRange, RectangularEtaPhiTrackingRegion::etaRange(), PixelRecoRange< T >::max(), PixelRecoRange< T >::mean(), CoreSimTrack::momentum(), PV3DBase< T, PVType, FrameType >::phi(), RectangularEtaPhiTrackingRegion::phiMargin(), TrackingRegionBase::ptMin(), PtMinSelector_cfg::ptMin, TkTrackingRegionsMargin< T >::right(), and mathSSE::sqrt().

Referenced by produce().

                                                      { 
    
  // 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());
  return inPtRange;
  
}

void FastTSGFromL2Muon::produce	(	edm::Event &	ev,
		const edm::EventSetup &	es
	)		`[virtual]`

Implements edm::EDProducer.

Definition at line 70 of file FastTSGFromL2Muon.cc.

References clean(), edm::Event::getByLabel(), iseed, edm::Event::put(), TrajectorySeed::recHits(), MuonTrackingRegionBuilder::region(), query::result, MuonTrackingRegionBuilder::setEvent(), SiTrackerGSMatchedRecHit2D::simtrackId(), theL2CollectionLabel, thePtCut, theRegionBuilder, theSeedCollectionLabels, theService, theSimTrackCollectionLabel, and MuonServiceProxy::update().

{

  // Initialize the output product
  std::auto_ptr<L3MuonTrajectorySeedCollection> result(new L3MuonTrajectorySeedCollection());

  //intialize service
  theService->update(es);
  
  // Region builder
  theRegionBuilder->setEvent(ev);

  // Retrieve the Monte Carlo truth (SimTracks)
  edm::Handle<edm::SimTrackContainer> theSimTracks;
  ev.getByLabel(theSimTrackCollectionLabel,theSimTracks);

  // Retrieve L2 muon collection
  edm::Handle<reco::TrackCollection> l2muonH;
  ev.getByLabel(theL2CollectionLabel ,l2muonH); 

  // Retrieve Seed collection
  unsigned seedCollections = theSeedCollectionLabels.size();
  std::vector<edm::Handle<edm::View<TrajectorySeed> > > theSeeds;
  theSeeds.resize(seedCollections);
  unsigned seed_size = 0;
  for ( unsigned iseed=0; iseed<seedCollections; ++iseed ) { 
    ev.getByLabel(theSeedCollectionLabels[iseed], theSeeds[iseed]);
    seed_size += theSeeds[iseed]->size();
  }

  // Loop on L2 muons
  unsigned int imu=0;
  unsigned int imuMax=l2muonH->size();
  // std::cout << "Found " << imuMax << " L2 muons" << std::endl;
  for (;imu!=imuMax;++imu){

    // Make a ref to l2 muon
    reco::TrackRef muRef(l2muonH, imu);
    
    // Cut on muons with low momenta
    if ( muRef->pt() < thePtCut 
         || muRef->innerMomentum().Rho() < thePtCut 
         || muRef->innerMomentum().R() < 2.5 ) continue;
    
    // Define the region of interest
    RectangularEtaPhiTrackingRegion * region = theRegionBuilder->region(muRef);

    // Copy the collection of seeds (ahem, this is time consuming!)
    std::vector<TrajectorySeed> tkSeeds;
    std::set<unsigned> tkIds;
    tkSeeds.reserve(seed_size);
    for ( unsigned iseed=0; iseed<seedCollections; ++iseed ) { 
      edm::Handle<edm::View<TrajectorySeed> > aSeedCollection = theSeeds[iseed];
      unsigned nSeeds = aSeedCollection->size();
      for (unsigned seednr = 0; seednr < nSeeds; ++seednr) {
            
        // The seed
        const BasicTrajectorySeed* aSeed = &((*aSeedCollection)[seednr]);
        
        // Find the first hit of the Seed
        TrajectorySeed::range theSeedingRecHitRange = aSeed->recHits();
        const SiTrackerGSMatchedRecHit2D * theFirstSeedingRecHit = 
          (const SiTrackerGSMatchedRecHit2D*) (&(*(theSeedingRecHitRange.first)));

        // The SimTrack id associated to that recHit
        int simTrackId = theFirstSeedingRecHit->simtrackId();

        // Track already associated to a seed
        std::set<unsigned>::iterator tkId = tkIds.find(simTrackId);
        if( tkId != tkIds.end() ) continue;     

        const SimTrack& theSimTrack = (*theSimTracks)[simTrackId]; 

        if ( clean(muRef,region,aSeed,theSimTrack) ) tkSeeds.push_back(*aSeed);
        tkIds.insert(simTrackId);

      } // End loop on seeds
 
    } // End loop on seed collections
  
    // Free memory
    delete region;

    // A plot
    // if(h_nSeedPerTrack) h_nSeedPerTrack->Fill(tkSeeds.size());

    // Another plot
    // if(h_nGoodSeedPerTrack) h_nGoodSeedPerTrack->Fill(tkSeeds.size());
    

    // Now create the Muon Trajectory Seed
    unsigned int is=0;
    unsigned int isMax=tkSeeds.size();
    for (;is!=isMax;++is){
      result->push_back( L3MuonTrajectorySeed(tkSeeds[is], muRef));
    } // End of tk seed loop
    
  } // End of l2 muon loop

  // std::cout << "Found " << result->size() << " seeds for muons" << std::endl;

  // And yet another plot
  // if(h_nGoodSeedPerEvent) h_nGoodSeedPerEvent->Fill(result->size());
  
  //put in the event
  ev.put(result);

}