MuonSeedCleaner.cc

Go to the documentation of this file.

#include <RecoMuon/MuonSeedGenerator/src/MuonSeedCleaner.h>

// Data Formats 
#include <DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h>
#include <DataFormats/TrajectorySeed/interface/TrajectorySeed.h>
#include <DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h>
#include <DataFormats/CSCRecHit/interface/CSCRecHit2D.h>
#include <DataFormats/CSCRecHit/interface/CSCSegmentCollection.h>
#include <DataFormats/CSCRecHit/interface/CSCSegment.h>
#include <DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h>
#include <DataFormats/DTRecHit/interface/DTRecSegment4D.h>

// Geometry
#include <Geometry/CommonDetUnit/interface/GeomDet.h>
#include <TrackingTools/DetLayers/interface/DetLayer.h>
#include <RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h>
#include <RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h>
#include <RecoMuon/Records/interface/MuonRecoGeometryRecord.h>

// muon service
#include <TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h>
#include <DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h>

// Framework
#include <FWCore/ParameterSet/interface/ParameterSet.h>
#include "FWCore/Framework/interface/Event.h"
#include <FWCore/Framework/interface/ESHandle.h>
#include <FWCore/MessageLogger/interface/MessageLogger.h> 
#include <DataFormats/Common/interface/Handle.h>

// C++
#include <vector>
#include <deque>
#include <utility>

//typedef std::pair<double, TrajectorySeed> seedpr ;
//static bool ptDecreasing(const seedpr s1, const seedpr s2) { return ( s1.first > s2.first ); }
static bool lengthSorting(const TrajectorySeed& s1, const TrajectorySeed& s2) { return ( s1.nHits() > s2.nHits() ); }

/*
 * Constructor
 */
MuonSeedCleaner::MuonSeedCleaner(const edm::ParameterSet& pset){

  // Local Debug flag
  debug                = pset.getParameter<bool>("DebugMuonSeed");

  // muon service
  edm::ParameterSet serviceParameters = pset.getParameter<edm::ParameterSet>("ServiceParameters");
  theService        = new MuonServiceProxy(serviceParameters);

}

/*
 * Destructor
 */
MuonSeedCleaner::~MuonSeedCleaner(){

  if (theService) delete theService;
}

/*********************************** 
 *
 *  Seed Cleaner
 *
 ***********************************/

std::vector<TrajectorySeed> MuonSeedCleaner::seedCleaner(const edm::EventSetup& eventSetup, std::vector<TrajectorySeed>& seeds ) {

  theService->update(eventSetup);
  
  std::vector<TrajectorySeed> FinalSeeds;
  
  // group the seeds
  std::vector<SeedContainer> theCollection = GroupSeeds(seeds);

  // ckeck each group and pick the good one
  for (size_t i=0; i< theCollection.size(); i++ ) {

      // separate seeds w/ more than 1 segments and w/ 1st layer segment information
      SeedContainer goodSeeds   = SeedCandidates( theCollection[i], true );
      SeedContainer otherSeeds  = SeedCandidates( theCollection[i], false );
      if ( MomentumFilter( goodSeeds ) )  {
         //std::cout<<" == type1 "<<std::endl;
         TrajectorySeed bestSeed = Chi2LengthSelection( goodSeeds );
     FinalSeeds.push_back( bestSeed );
           
         GlobalPoint seedgp = SeedPosition( bestSeed );
         double eta = fabs( seedgp.eta() );
         if ( goodSeeds.size() > 2 && eta > 1.5 ) {
            TrajectorySeed anotherSeed  = MoreRecHits( goodSeeds );
            FinalSeeds.push_back( anotherSeed );
         }
      } 
      else if( MomentumFilter( otherSeeds ) ) {
         //std::cout<<" == type2 "<<std::endl;
         TrajectorySeed bestSeed = MoreRecHits( otherSeeds );
     FinalSeeds.push_back( bestSeed );

         GlobalPoint seedgp = SeedPosition( bestSeed );
         double eta = fabs( seedgp.eta() );
         if ( otherSeeds.size() > 2 && eta > 1.5 ) {
            TrajectorySeed anotherSeed  = LeanHighMomentum( otherSeeds );
            FinalSeeds.push_back( anotherSeed );
         }
      } 
      else {
         //std::cout<<" == type3 "<<std::endl;
         TrajectorySeed bestSeed = LeanHighMomentum( theCollection[i] );
     FinalSeeds.push_back( bestSeed );

         GlobalPoint seedgp = SeedPosition( bestSeed );
         double eta = fabs( seedgp.eta() );
         if ( theCollection.size() > 2 && eta > 1.5 ) {
            TrajectorySeed anotherSeed  = BiggerCone( theCollection[i] );
            FinalSeeds.push_back( anotherSeed );
         }
      }
  }  
  return FinalSeeds ; 

}


TrajectorySeed MuonSeedCleaner::Chi2LengthSelection(std::vector<TrajectorySeed>& seeds ) {

  if ( seeds.size() == 1 ) return seeds[0];

  int    winner   = 0 ;
  int    moreHits = 0 ;  
  double bestChi2 = 99999.;  
  for ( size_t i = 0; i < seeds.size(); i++ ) {

      // 1. fill out the Nchi2 of segments of the seed 
      //GlobalVector mom = SeedMomentum( seeds[i] ); // temporary use for debugging
      //double pt = sqrt( (mom.x()*mom.x()) + (mom.y()*mom.y()) );
      //std::cout<<" > SEED"<<i<<"  pt:"<<pt<< std::endl;

      double theChi2 = SeedChi2( seeds[i] );  
      double dChi2 = fabs( 1. - (theChi2 / bestChi2)) ;
      int    theHits = seeds[i].nHits() ;  
      int    dHits = theHits - moreHits ;
      //std::cout<<" -----  "<<std::endl;

      // 2. better chi2 
      if ( theChi2 < bestChi2 && dChi2 > 0.05 ) {
           winner = static_cast<int>(i) ;
           bestChi2 = theChi2 ;
           moreHits = theHits ;
      }
      // 3. if chi2 is not much better, pick more rechits one 
      if ( theChi2 >= bestChi2 && dChi2 < 0.05 && dHits > 0 ) {
           winner = static_cast<int>(i) ;
           bestChi2 = theChi2 ;
           moreHits = theHits ;
      }

  }
  //std::cout<<" Winner is "<< winner <<std::endl;
  TrajectorySeed theSeed =  seeds[winner];
  seeds.erase( seeds.begin()+winner ); 
  return theSeed;
}

TrajectorySeed MuonSeedCleaner::BiggerCone(std::vector<TrajectorySeed>& seeds ) {

  if ( seeds.size() == 1 ) return seeds[0];

  float biggerProjErr  = 9999.; 
  int winner = 0 ;
  AlgebraicSymMatrix mat(5,0) ; 
  for ( size_t i = 0; i < seeds.size(); i++ ) {

      edm::OwnVector<TrackingRecHit>::const_iterator r1 = seeds[i].recHits().first ;
      mat = r1->parametersError().similarityT( r1->projectionMatrix() );

      int    NRecHits = NRecHitsFromSegment( *r1 );

      float ddx = mat[1][1]; 
      float ddy = mat[2][2]; 
      float dxx = mat[3][3]; 
      float dyy = mat[4][4];
      float projectErr = sqrt( (ddx*10000.) + (ddy*10000.) + dxx + dyy ) ;

      if ( NRecHits < 5 ) continue ;
      if ( projectErr < biggerProjErr ) continue;

      winner = static_cast<int>(i) ;
      biggerProjErr = projectErr ;
  }
  TrajectorySeed theSeed = seeds[winner];
  seeds.erase( seeds.begin()+winner ); 
  return theSeed;

}

TrajectorySeed MuonSeedCleaner::LeanHighMomentum( std::vector<TrajectorySeed>& seeds ) {

  if ( seeds.size() == 1 ) return seeds[0];

  double highestPt = 0. ;  
  int    winner    = 0  ;
  for ( size_t i = 0; i < seeds.size(); i++ ) {
       GlobalVector mom = SeedMomentum( seeds[i] );
       double pt = sqrt( (mom.x()*mom.x()) + (mom.y()*mom.y()) );
       if ( pt > highestPt ) { 
          winner = static_cast<int>(i) ;
          highestPt = pt ;
       }
  }
  TrajectorySeed theSeed = seeds[winner];
  seeds.erase( seeds.begin()+winner ); 
  return theSeed;

}


TrajectorySeed MuonSeedCleaner::MoreRecHits(std::vector<TrajectorySeed>& seeds ) {

  if ( seeds.size() == 1 ) return seeds[0];

  int    winner   = 0 ;
  int    moreHits = 0 ;  
  double betterChi2 = 99999.;
  for ( size_t i = 0; i < seeds.size(); i++ ) {

      int    theHits = 0;  
      for (edm::OwnVector<TrackingRecHit>::const_iterator r1 = seeds[i].recHits().first; r1 != seeds[i].recHits().second; r1++){
          theHits += NRecHitsFromSegment( *r1 );
      }

      double theChi2 = SeedChi2( seeds[i] );  
   
      if ( theHits == moreHits && theChi2 < betterChi2  ) {
         betterChi2 = theChi2; 
         winner = static_cast<int>(i) ;
      } 
      if ( theHits > moreHits ) {
           moreHits = theHits ;
           betterChi2 = theChi2;
           winner = static_cast<int>(i) ;
      }
  }
  TrajectorySeed theSeed = seeds[winner];
  seeds.erase( seeds.begin()+winner ); 
  return theSeed;
}


SeedContainer MuonSeedCleaner::LengthFilter(std::vector<TrajectorySeed>& seeds ) {
 
  SeedContainer longSeeds; 
  int NSegs = 0;
  for (size_t i = 0; i< seeds.size(); i++) {
    
      int theLength = static_cast<int>( seeds[i].nHits());
      if ( theLength > NSegs ) {
         NSegs = theLength ;
         longSeeds.clear();
         longSeeds.push_back( seeds[i] );
      } 
      else if ( theLength == NSegs ) {
         longSeeds.push_back( seeds[i] );
      } else {
         continue;
      } 
  }
  //std::cout<<" final Length :"<<NSegs<<std::endl;

  return longSeeds ; 
  
}


bool MuonSeedCleaner::MomentumFilter(std::vector<TrajectorySeed>& seeds ) {

  bool findgoodMomentum = false;
  SeedContainer goodMomentumSeeds = seeds;
  seeds.clear();
  for ( size_t i = 0; i < goodMomentumSeeds.size(); i++ ) {
       GlobalVector mom = SeedMomentum( goodMomentumSeeds[i] );
       double pt = sqrt( (mom.x()*mom.x()) + (mom.y()*mom.y()) );
       if ( pt < 6. || pt > 2000. ) continue;
       //if ( pt < 6. ) continue;
       //std::cout<<" passed momentum :"<< pt <<std::endl;
       seeds.push_back( goodMomentumSeeds[i] );
       findgoodMomentum = true;  
  }
  if ( seeds.empty() ) seeds = goodMomentumSeeds;

  return findgoodMomentum;
}



SeedContainer MuonSeedCleaner::SeedCandidates( std::vector<TrajectorySeed>& seeds, bool good ) {

  SeedContainer theCandidate;
  theCandidate.clear();

  bool longSeed = false;  
  bool withFirstLayer = false ;

  //std::cout<<"***** Seed Classification *****"<< seeds.size() <<std::endl;
  for ( size_t i = 0; i < seeds.size(); i++ ) {

      if (seeds[i].nHits() > 1 ) longSeed = true ;
      //std::cout<<"  Seed: "<<i<<" w/"<<seeds[i].nHits()<<" segs "<<std::endl;
      // looking for 1st layer segment
      int idx = 0;
      for (edm::OwnVector<TrackingRecHit>::const_iterator r1 = seeds[i].recHits().first; r1 != seeds[i].recHits().second; r1++){

         idx++;
         const GeomDet* gdet = theService->trackingGeometry()->idToDet( (*r1).geographicalId() );
         DetId geoId = gdet->geographicalId();

         if ( geoId.subdetId() == MuonSubdetId::DT ) {
            DTChamberId DT_Id( (*r1).geographicalId() );
            //std::cout<<" ID:"<<DT_Id <<" pos:"<< r1->localPosition()  <<std::endl;
            if (DT_Id.station() != 1)  continue;
            withFirstLayer = true;
         }
         if ( geoId.subdetId() == MuonSubdetId::CSC ) {
            idx++;
            CSCDetId CSC_Id = CSCDetId( (*r1).geographicalId() );
            //std::cout<<" ID:"<<CSC_Id <<" pos:"<< r1->localPosition()  <<std::endl;
            if (CSC_Id.station() != 1)  continue;
            withFirstLayer = true;
         }
      }
      bool goodseed = (longSeed && withFirstLayer) ? true : false ;
  
      if ( goodseed == good )  theCandidate.push_back( seeds[i] );
  }
  return theCandidate;

}

std::vector<SeedContainer> MuonSeedCleaner::GroupSeeds( std::vector<TrajectorySeed>& seeds) {

  std::vector<SeedContainer> seedCollection;
  seedCollection.clear();
  std::vector<TrajectorySeed> theGroup ;
  std::vector<bool> usedSeed(seeds.size(),false);

  // categorize seeds by comparing overlapping segments or a certian eta-phi cone 
  for (size_t i= 0; i<seeds.size(); i++){
    
    if (usedSeed[i]) continue;
    theGroup.push_back( seeds[i] );
    usedSeed[i] = true ;

    GlobalPoint pos1 = SeedPosition( seeds[i] );

    for (size_t j= i+1; j<seeds.size(); j++){
 
       // 1.1 seeds with overlaaping segments will be grouped together
       unsigned int overlapping = OverlapSegments(seeds[i], seeds[j]) ;
       if ( !usedSeed[j] && overlapping > 0 ) {
          // reject the identical seeds
          if ( seeds[i].nHits() == overlapping && seeds[j].nHits() == overlapping ) {
             usedSeed[j] = true ;
             continue;
          }
          theGroup.push_back( seeds[j] ); 
          usedSeed[j] = true ;
       }
       if (usedSeed[j]) continue;

       // 1.2 seeds in a certain cone are grouped together  
       GlobalPoint pos2 = SeedPosition( seeds[j]);
       double dh = pos1.eta() - pos2.eta() ;
       double df = pos1.phi() - pos2.phi() ;
       double dR = sqrt( (dh*dh) + (df*df) );

       if ( dR > 0.3 && seeds[j].nHits() == 1) continue;
       if ( dR > 0.2 && seeds[j].nHits() >  1) continue;
       theGroup.push_back( seeds[j] ); 
       usedSeed[j] = true ;
    }
    sort(theGroup.begin(), theGroup.end(), lengthSorting ) ;
    seedCollection.push_back(theGroup);
    //std::cout<<" group "<<seedCollection.size() <<" w/"<< theGroup.size() <<" seeds"<<std::endl; 
    theGroup.clear(); 
  }
  return seedCollection;

}

unsigned int MuonSeedCleaner::OverlapSegments( const TrajectorySeed& seed1, const TrajectorySeed& seed2 ) {

  unsigned int overlapping = 0;
  for (edm::OwnVector<TrackingRecHit>::const_iterator r1 = seed1.recHits().first; r1 != seed1.recHits().second; r1++){
      
      DetId id1 = (*r1).geographicalId();
      const GeomDet* gdet1 = theService->trackingGeometry()->idToDet( id1 );
      GlobalPoint gp1 = gdet1->toGlobal( (*r1).localPosition() );

      for (edm::OwnVector<TrackingRecHit>::const_iterator r2 = seed2.recHits().first; r2 != seed2.recHits().second; r2++){

          DetId id2 = (*r2).geographicalId();
          if (id1 != id2 ) continue;

      const GeomDet* gdet2 = theService->trackingGeometry()->idToDet( id2 );
      GlobalPoint gp2 = gdet2->toGlobal( (*r2).localPosition() );

      double dx = gp1.x() - gp2.x() ;
      double dy = gp1.y() - gp2.y() ;
      double dz = gp1.z() - gp2.z() ;
      double dL = sqrt( dx*dx + dy*dy + dz*dz);

          if ( dL < 1. ) overlapping ++;
      
      }
  }
  return overlapping ;

}

double MuonSeedCleaner::SeedChi2( const TrajectorySeed& seed ) {

  double theChi2 = 0.;  
  for (edm::OwnVector<TrackingRecHit>::const_iterator r1 = seed.recHits().first; r1 != seed.recHits().second; r1++){
      //std::cout<<"    segmet : "<<it <<std::endl; 
      theChi2 += NChi2OfSegment( *r1 );
  }
  theChi2 = theChi2 / seed.nHits() ;

  //std::cout<<" final Length :"<<NSegs<<std::endl;
  return theChi2 ;

}

int MuonSeedCleaner::SeedLength( const TrajectorySeed& seed ) {

  int theHits = 0;  
  for (edm::OwnVector<TrackingRecHit>::const_iterator r1 = seed.recHits().first; r1 != seed.recHits().second; r1++){
      //std::cout<<"    segmet : "<<it <<std::endl; 
      theHits += NRecHitsFromSegment( *r1 );
  }

  //std::cout<<" final Length :"<<NSegs<<std::endl;
  return theHits ;

}

GlobalPoint MuonSeedCleaner::SeedPosition( const TrajectorySeed& seed ) {

  

  PTrajectoryStateOnDet pTSOD = seed.startingState();
  DetId SeedDetId(pTSOD.detId());
  const GeomDet* geoDet = theService->trackingGeometry()->idToDet( SeedDetId );
  TrajectoryStateOnSurface SeedTSOS = trajectoryStateTransform::transientState(pTSOD, &(geoDet->surface()), &*theService->magneticField());
  GlobalPoint  pos  = SeedTSOS.globalPosition();

  return pos ;

}

GlobalVector MuonSeedCleaner::SeedMomentum( const TrajectorySeed& seed ) {

  

  PTrajectoryStateOnDet pTSOD = seed.startingState();
  DetId SeedDetId(pTSOD.detId());
  const GeomDet* geoDet = theService->trackingGeometry()->idToDet( SeedDetId );
  TrajectoryStateOnSurface SeedTSOS = trajectoryStateTransform::transientState(pTSOD, &(geoDet->surface()), &*theService->magneticField());
  GlobalVector  mom  = SeedTSOS.globalMomentum();

  return mom ;

}

int MuonSeedCleaner::NRecHitsFromSegment( const TrackingRecHit& rhit ) {

      int NRechits = 0 ; 
      const GeomDet* gdet = theService->trackingGeometry()->idToDet( rhit.geographicalId() );
      MuonTransientTrackingRecHit::MuonRecHitPointer theSeg = 
      MuonTransientTrackingRecHit::specificBuild(gdet, rhit.clone() );

      DetId geoId = gdet->geographicalId();
      if ( geoId.subdetId() == MuonSubdetId::DT ) {
         DTChamberId DT_Id( rhit.geographicalId() );
     std::vector<TrackingRecHit*> DThits = theSeg->recHits();
     int dt1DHits = 0;
     for (size_t j=0; j< DThits.size(); j++) {
             dt1DHits += (DThits[j]->recHits()).size();
         }
         NRechits = dt1DHits ;
      }

      if ( geoId.subdetId() == MuonSubdetId::CSC ) {
         NRechits = (theSeg->recHits()).size() ;
      }
      return NRechits ;
}

int MuonSeedCleaner::NRecHitsFromSegment( MuonTransientTrackingRecHit *rhit ) {

    int NRechits = 0 ; 
    DetId geoId = rhit->geographicalId();
    if ( geoId.subdetId() == MuonSubdetId::DT ) {
       DTChamberId DT_Id( geoId );
       std::vector<TrackingRecHit*> DThits = rhit->recHits();
       int dt1DHits = 0;
       for (size_t j=0; j< DThits.size(); j++) {
           dt1DHits += (DThits[j]->recHits()).size();
       }
       NRechits = dt1DHits ;
       //std::cout<<" D_rh("<< dt1DHits  <<") " ;
    }
    if ( geoId.subdetId() == MuonSubdetId::CSC ) {
       NRechits = (rhit->recHits()).size() ;
       //std::cout<<" C_rh("<<(rhit->recHits()).size() <<") " ;
    }
    return NRechits;

}

double MuonSeedCleaner::NChi2OfSegment( const TrackingRecHit& rhit ) {

      double NChi2 = 999999. ; 
      const GeomDet* gdet = theService->trackingGeometry()->idToDet( rhit.geographicalId() );
      MuonTransientTrackingRecHit::MuonRecHitPointer theSeg = 
      MuonTransientTrackingRecHit::specificBuild(gdet, rhit.clone() );

      double dof = static_cast<double>( theSeg->degreesOfFreedom() );
      NChi2 = theSeg->chi2() / dof ;
      //std::cout<<" Chi2 = "<< NChi2  <<" |" ;

      return NChi2 ;
}