TrackTransformerForGlobalCosmicMuons.cc

Go to the documentation of this file.

#include "TrackingTools/TrackRefitter/interface/TrackTransformerForGlobalCosmicMuons.h"

#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"

#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
#include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
#include "TrackingTools/PatternTools/interface/TrajectorySmoother.h"

#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
#include "TrackingTools/Records/interface/TransientRecHitRecord.h"
#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "TrackingTools/GeomPropagators/interface/Propagator.h"

#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/DetId/interface/DetId.h"


#include "DataFormats/SiStripDetId/interface/TIBDetId.h"
#include "DataFormats/SiStripDetId/interface/TOBDetId.h"
#include "DataFormats/SiStripDetId/interface/TIDDetId.h"
#include "DataFormats/SiStripDetId/interface/TECDetId.h"
#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
#include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
#include "DataFormats/MuonDetId/interface/CSCDetId.h"
#include "DataFormats/MuonDetId/interface/DTWireId.h"
#include "DataFormats/MuonDetId/interface/RPCDetId.h"
#include "DataFormats/MuonDetId/interface/MuonSubdetId.h"


using namespace std;
using namespace edm;

TrackTransformerForGlobalCosmicMuons::TrackTransformerForGlobalCosmicMuons(const ParameterSet& parameterSet){
  
  theTrackerRecHitBuilderName = parameterSet.getParameter<string>("TrackerRecHitBuilder");
  theMuonRecHitBuilderName = parameterSet.getParameter<string>("MuonRecHitBuilder");

  theRPCInTheFit = parameterSet.getParameter<bool>("RefitRPCHits");

  theCacheId_TC = theCacheId_GTG = theCacheId_MG = theCacheId_TRH = 0;
  theSkipStationDT      = parameterSet.getParameter<int>("SkipStationDT");
  theSkipStationCSC     = parameterSet.getParameter<int>("SkipStationCSC");
  theSkipWheelDT        = parameterSet.getParameter<int>("SkipWheelDT");
  theTrackerSkipSystem  = parameterSet.getParameter<int>("TrackerSkipSystem");
  theTrackerSkipSection = parameterSet.getParameter<int>("TrackerSkipSection");//layer, wheel, or disk depending on the system 
}

TrackTransformerForGlobalCosmicMuons::~TrackTransformerForGlobalCosmicMuons(){}


void TrackTransformerForGlobalCosmicMuons::setServices(const EventSetup& setup){
  
  const std::string metname = "Reco|TrackingTools|TrackTransformer";

  setup.get<TrajectoryFitter::Record>().get("KFFitterForRefitInsideOut",theFitterIO);
  setup.get<TrajectoryFitter::Record>().get("KFSmootherForRefitInsideOut",theSmootherIO);
  setup.get<TrajectoryFitter::Record>().get("KFFitterForRefitOutsideIn",theFitterOI);
  setup.get<TrajectoryFitter::Record>().get("KFSmootherForRefitOutsideIn",theSmootherOI);
  
  unsigned long long newCacheId_TC = setup.get<TrackingComponentsRecord>().cacheIdentifier();

  if ( newCacheId_TC != theCacheId_TC ){
    LogTrace(metname) << "Tracking Component changed!";
    theCacheId_TC = newCacheId_TC;
    setup.get<TrackingComponentsRecord>().get("SmartPropagatorRK",thePropagatorIO);
    setup.get<TrackingComponentsRecord>().get("SmartPropagatorRKOpposite",thePropagatorOI);

  }

  // Global Tracking Geometry
  unsigned long long newCacheId_GTG = setup.get<GlobalTrackingGeometryRecord>().cacheIdentifier();
  if ( newCacheId_GTG != theCacheId_GTG ) {
    LogTrace(metname) << "GlobalTrackingGeometry changed!";
    theCacheId_GTG = newCacheId_GTG;
    setup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry); 
  }
  
  // Magfield Field
  unsigned long long newCacheId_MG = setup.get<IdealMagneticFieldRecord>().cacheIdentifier();
  if ( newCacheId_MG != theCacheId_MG ) {
    LogTrace(metname) << "Magnetic Field changed!";
    theCacheId_MG = newCacheId_MG;
    setup.get<IdealMagneticFieldRecord>().get(theMGField);
  }
  
  // Transient Rechit Builders
  unsigned long long newCacheId_TRH = setup.get<TransientRecHitRecord>().cacheIdentifier();
  if ( newCacheId_TRH != theCacheId_TRH ) {
    theCacheId_TRH = newCacheId_TRH;
    LogTrace(metname) << "TransientRecHitRecord changed!";
    setup.get<TransientRecHitRecord>().get(theTrackerRecHitBuilderName,theTrackerRecHitBuilder);
    setup.get<TransientRecHitRecord>().get(theMuonRecHitBuilderName,theMuonRecHitBuilder);
  }
}


TransientTrackingRecHit::ConstRecHitContainer
TrackTransformerForGlobalCosmicMuons::getTransientRecHits(const reco::TransientTrack& track) const {

  TransientTrackingRecHit::ConstRecHitContainer tkHits;
  TransientTrackingRecHit::ConstRecHitContainer staHits;

  for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
    if((*hit)->isValid()) {
      if ( (*hit)->geographicalId().det() == DetId::Tracker && TrackerKeep((*hit)->geographicalId())) {
        tkHits.push_back(theTrackerRecHitBuilder->build(&**hit));
      } else if ( (*hit)->geographicalId().det() == DetId::Muon && MuonKeep((*hit)->geographicalId())){
        if( (*hit)->geographicalId().subdetId() == 3 && !theRPCInTheFit){
        LogTrace("Reco|TrackingTools|TrackTransformer") << "RPC Rec Hit discarged"; 
        continue;
        }
        staHits.push_back(theMuonRecHitBuilder->build(&**hit));
      }
    }
  }
  
  if(staHits.empty()) return staHits;

/*
  bool up = staHits.front()->globalPosition().y()>0 ? true : false;

  if(up){
    reverse(staHits.begin(),staHits.end());
    reverse(tkHits.begin(),tkHits.end());
  }
*/
  copy(staHits.begin(),staHits.end(),back_inserter(tkHits));

  for(TransientTrackingRecHit::ConstRecHitContainer::const_iterator hit = tkHits.begin();
      hit !=tkHits.end(); ++hit){

    DetId hitId = (*hit)->geographicalId();
    GlobalPoint glbpoint = trackingGeometry()->idToDet(hitId)->position();

    if(hitId.det() == DetId::Tracker) {
      if (hitId.subdetId() == StripSubdetector::TIB )  
    LogTrace("TrackFitters") << glbpoint << " I am TIB " << TIBDetId(hitId).layer();
      else if (hitId.subdetId() == StripSubdetector::TOB ) 
    LogTrace("TrackFitters") << glbpoint << " I am TOB " << TOBDetId(hitId).layer();
      else if (hitId.subdetId() == StripSubdetector::TEC ) 
    LogTrace("TrackFitters") << glbpoint << " I am TEC " << TECDetId(hitId).wheel();
      else if (hitId.subdetId() == StripSubdetector::TID ) 
    LogTrace("TrackFitters") << glbpoint << " I am TID " << TIDDetId(hitId).wheel();
      else if (hitId.subdetId() == StripSubdetector::TID ) 
    LogTrace("TrackFitters") << glbpoint << " I am TID " << TIDDetId(hitId).wheel();
      else if (hitId.subdetId() == (int) PixelSubdetector::PixelBarrel ) 
    LogTrace("TrackFitters") << glbpoint << " I am PixBar " << PXBDetId(hitId).layer();
      else if (hitId.subdetId() == (int) PixelSubdetector::PixelEndcap )
    LogTrace("TrackFitters") << glbpoint << " I am PixFwd " << PXFDetId(hitId).disk();
      else 
    LogTrace("TrackFitters") << " UNKNOWN TRACKER HIT TYPE ";
    } else if(hitId.det() == DetId::Muon) {
      if(hitId.subdetId() == MuonSubdetId::DT)
    LogTrace("TrackFitters") << glbpoint << " I am DT " << DTWireId(hitId);
      else if (hitId.subdetId() == MuonSubdetId::CSC )
    LogTrace("TrackFitters") << glbpoint << " I am CSC " << CSCDetId(hitId);
      else if (hitId.subdetId() == MuonSubdetId::RPC )
    LogTrace("TrackFitters") << glbpoint << " I am RPC " << RPCDetId(hitId);
      else 
    LogTrace("TrackFitters") << " UNKNOWN MUON HIT TYPE ";
    } else
      LogTrace("TrackFitters") << " UNKNOWN HIT TYPE ";
  } 
  
  return tkHits;
}


ESHandle<TrajectoryFitter> TrackTransformerForGlobalCosmicMuons::fitter(bool up) const{
  if(up) return theFitterOI;
  else return theFitterIO;
}
  
ESHandle<TrajectorySmoother> TrackTransformerForGlobalCosmicMuons::smoother(bool up) const{
  if(up) return theSmootherOI;
  else return theSmootherIO;
}

ESHandle<Propagator> TrackTransformerForGlobalCosmicMuons::propagator(bool up) const{
  if(up) return thePropagatorIO;
  else return thePropagatorOI;
}



vector<Trajectory> TrackTransformerForGlobalCosmicMuons::transform(const reco::Track& tr) const {

  const std::string metname = "Reco|TrackingTools|TrackTransformer";
  
  reco::TransientTrack track(tr,magneticField(),trackingGeometry());   

  // Build the transient Rechits
  TransientTrackingRecHit::ConstRecHitContainer recHitsForReFit = getTransientRecHits(track);

  if(recHitsForReFit.size() < 2) return vector<Trajectory>();

  bool up = recHitsForReFit.back()->globalPosition().y()>0 ? true : false;
  LogTrace(metname) << "Up ? " << up;

  PropagationDirection propagationDirection = up ? oppositeToMomentum : alongMomentum;
  TrajectoryStateOnSurface firstTSOS = up ? track.outermostMeasurementState() : track.innermostMeasurementState();
  unsigned int innerId = up ? track.track().outerDetId() : track.track().innerDetId();

  LogTrace(metname) << "Prop Dir: " << propagationDirection << " FirstId " << innerId << " firstTSOS " << firstTSOS;

  TrajectorySeed seed(PTrajectoryStateOnDet(),TrajectorySeed::recHitContainer(),propagationDirection);


  if(recHitsForReFit.front()->geographicalId() != DetId(innerId)){
    LogTrace(metname)<<"Propagation occurring"<<endl;
    firstTSOS = propagator(up)->propagate(firstTSOS, recHitsForReFit.front()->det()->surface());
    LogTrace(metname)<<"Final destination: " << recHitsForReFit.front()->det()->surface().position() << endl;
    if(!firstTSOS.isValid()){
      LogTrace(metname)<<"Propagation error!"<<endl;
      return vector<Trajectory>();
    }
  }
  

  vector<Trajectory> trajectories = fitter(up)->fit(seed,recHitsForReFit,firstTSOS);
  
  if(trajectories.empty()){
    LogTrace(metname)<<"No Track refitted!"<<endl;
    return vector<Trajectory>();
  }
  
  Trajectory trajectoryBW = trajectories.front();
    
  vector<Trajectory> trajectoriesSM = smoother(up)->trajectories(trajectoryBW);

  if(trajectoriesSM.empty()){
    LogTrace(metname)<<"No Track smoothed!"<<endl;
    return vector<Trajectory>();
  }
  
  return trajectoriesSM;

}

//
// Selection for Tracker Hits
//
bool TrackTransformerForGlobalCosmicMuons::TrackerKeep(DetId id) const{
    
    bool retVal = true;
    if (id.det() != DetId::Tracker ) return false;
    if (theTrackerSkipSystem < 0 ) return true;
    

    int layer = -999;
    int disk  = -999;
    int wheel = -999;

      if ( id.subdetId() == theTrackerSkipSystem){

        if (theTrackerSkipSystem == PixelSubdetector::PixelBarrel) {
            PXBDetId did(id.rawId());
            layer = did.layer();
        }

        if (theTrackerSkipSystem == StripSubdetector::TIB) {
            TIBDetId did(id.rawId());
            layer = did.layer();
        }   

        if (theTrackerSkipSystem == StripSubdetector::TOB) {
            TOBDetId did(id.rawId());
            layer = did.layer();
        }
        if (theTrackerSkipSystem == PixelSubdetector::PixelEndcap) {
            PXFDetId did(id.rawId());
            disk = did.disk();
        }
        if (theTrackerSkipSystem == StripSubdetector::TID) {
            TIDDetId did(id.rawId());
            wheel = did.wheel();
        }
        if (theTrackerSkipSystem == StripSubdetector::TEC) {
            TECDetId did(id.rawId());
            wheel = did.wheel();
        }
    }

    if (theTrackerSkipSection > -998 && layer == theTrackerSkipSection) retVal = false;
    if (theTrackerSkipSection > -998 && disk  == theTrackerSkipSection) retVal = false;
    if (theTrackerSkipSection > -998 && wheel == theTrackerSkipSection) retVal = false;

    return retVal;
}
//
// Selection for Muon hits
//
bool TrackTransformerForGlobalCosmicMuons::MuonKeep(DetId id) const {

    if (id.det() != DetId::Muon) return false;
    if (theSkipStationDT < 0 && theSkipStationCSC < 0) return true;

    int station = -999;
    int wheel   = -999;
    bool isRPC  = false;
    bool isDT   = false;
    bool isCSC  = false;
        
    if ( id.subdetId() == MuonSubdetId::DT ) {
        DTChamberId did(id.rawId());
        station = did.station();
        wheel = did.wheel();
        isDT = true;
    } else if ( id.subdetId() == MuonSubdetId::CSC ) {
        CSCDetId did(id.rawId());
        station = did.station();
        isCSC = true;
    } else if ( id.subdetId() == MuonSubdetId::RPC ) {
        RPCDetId rpcid(id.rawId());
        station = rpcid.station();
        isRPC = true;
    }
    
    if (isRPC   && (station == theSkipStationCSC || station == theSkipStationDT)) return false;
    if (isDT    && station  == theSkipStationDT     ) return false;
    if (isCSC   && station  == theSkipStationCSC    ) return false;

    if (isDT && theSkipWheelDT > -998 && wheel == theSkipWheelDT) return false;
    
    return true;
}