KFTrajectorySmoother Class Reference

#include <KFTrajectorySmoother.h>

Inheritance diagram for KFTrajectorySmoother:

Public Member Functions
const Propagator *	alongPropagator () const
virtual KFTrajectorySmoother *	clone () const override
const MeasurementEstimator *	estimator () const
	KFTrajectorySmoother (const Propagator &aPropagator, const TrajectoryStateUpdator &aUpdator, const MeasurementEstimator &aEstimator, float errorRescaling=100.f, int minHits=3)
	KFTrajectorySmoother (const Propagator *aPropagator, const TrajectoryStateUpdator *aUpdator, const MeasurementEstimator *aEstimator, float errorRescaling=100.f, int minHits=3, const DetLayerGeometry *detLayerGeometry=0, TkCloner const *hc=0)
const Propagator *	oppositePropagator () const
virtual void	setHitCloner (TkCloner const *hc)
virtual Trajectory	trajectory (const Trajectory &aTraj) const override
const TrajectoryStateUpdator *	updator () const
virtual	~KFTrajectorySmoother ()
Public Member Functions inherited from TrajectorySmoother
virtual TrajectoryContainer	trajectories (const Trajectory &traj) const
virtual	~TrajectorySmoother ()

Private Types
typedef FreeTrajectoryState	FTS
typedef TrajectoryMeasurement	TM
typedef TrajectoryStateOnSurface	TSOS

Private Attributes
const DetLayerGeometry	dummyGeometry
int	minHits_
const Propagator *	theAlongPropagator
float	theErrorRescaling
const MeasurementEstimator *	theEstimator
const DetLayerGeometry *	theGeometry
TkCloner const *	theHitCloner = 0
const Propagator *	theOppositePropagator
const TrajectoryStateUpdator *	theUpdator

Additional Inherited Members
Public Types inherited from TrajectorySmoother
typedef std::vector< Trajectory >	TrajectoryContainer
typedef TrajectoryContainer::iterator	TrajectoryIterator

Detailed Description

A Standard Kalman smoother. The forward fit is not redone, only the backward smoothing. Ported from ORCA

Author

todorov, cerati

Definition at line 20 of file KFTrajectorySmoother.h.

Member Typedef Documentation

typedef FreeTrajectoryState KFTrajectorySmoother::FTS

private

Definition at line 25 of file KFTrajectorySmoother.h.

typedef TrajectoryMeasurement KFTrajectorySmoother::TM

private

Definition at line 26 of file KFTrajectorySmoother.h.

typedef TrajectoryStateOnSurface KFTrajectorySmoother::TSOS

private

Definition at line 24 of file KFTrajectorySmoother.h.

Constructor & Destructor Documentation

KFTrajectorySmoother::KFTrajectorySmoother ( const Propagator &  aPropagator, const TrajectoryStateUpdator &  aUpdator, const MeasurementEstimator &  aEstimator, float  errorRescaling = 100.f, int  minHits = 3  )

inline

Definition at line 30 of file KFTrajectorySmoother.h.

References alongMomentum, Propagator::clone(), dummyGeometry, oppositeToMomentum, AlCaHLTBitMon_ParallelJobs::p, theAlongPropagator, theGeometry, and theOppositePropagator.

Referenced by clone().

                                :
  theAlongPropagator(nullptr),
    theOppositePropagator(nullptr),
    theUpdator(aUpdator.clone()),
    theEstimator(aEstimator.clone()),
    theErrorRescaling(errorRescaling),
    minHits_(minHits),
    theGeometry(nullptr){ // to be fixed. Why this first constructor is needed? who is using it? Can it be removed?
      if(!theGeometry) theGeometry = &dummyGeometry;
      auto p = aPropagator.clone();
      p->setPropagationDirection(alongMomentum);
      theAlongPropagator = p;
      p = aPropagator.clone();
      p->setPropagationDirection(oppositeToMomentum);
      theOppositePropagator = p;
    }

KFTrajectorySmoother::KFTrajectorySmoother ( const Propagator *  aPropagator, const TrajectoryStateUpdator *  aUpdator, const MeasurementEstimator *  aEstimator, float  errorRescaling = 100.f, int  minHits = 3, const DetLayerGeometry *  detLayerGeometry = 0, TkCloner const *  hc = 0  )

inline

Definition at line 52 of file KFTrajectorySmoother.h.

References alongMomentum, Propagator::clone(), dummyGeometry, oppositeToMomentum, AlCaHLTBitMon_ParallelJobs::p, theAlongPropagator, theGeometry, and theOppositePropagator.

                                                    :
    theAlongPropagator(nullptr),
    theOppositePropagator(nullptr),
    theUpdator(aUpdator->clone()),
    theEstimator(aEstimator->clone()),
    theHitCloner(hc),
    theErrorRescaling(errorRescaling),
    minHits_(minHits),
    theGeometry(detLayerGeometry){
      if(!theGeometry) theGeometry = &dummyGeometry;
      auto p = aPropagator->clone();
      p->setPropagationDirection(alongMomentum);
      theAlongPropagator = p;
      p = aPropagator->clone();
      p->setPropagationDirection(oppositeToMomentum);
      theOppositePropagator = p;
    }

KFTrajectorySmoother::~KFTrajectorySmoother ( )

virtual

Definition at line 26 of file KFTrajectorySmoother.cc.

References theAlongPropagator, theEstimator, theOppositePropagator, and theUpdator.

                                            {

  delete theAlongPropagator;
  delete theOppositePropagator;
  delete theUpdator;
  delete theEstimator;

}

Member Function Documentation

const Propagator* KFTrajectorySmoother::alongPropagator ( ) const

inline

Definition at line 80 of file KFTrajectorySmoother.h.

References theAlongPropagator.

{ return theAlongPropagator;}

virtual KFTrajectorySmoother* KFTrajectorySmoother::clone ( void  ) const

inlineoverridevirtual

Implements TrajectorySmoother.

Definition at line 86 of file KFTrajectorySmoother.h.

References KFTrajectorySmoother(), minHits_, theAlongPropagator, theErrorRescaling, theEstimator, theGeometry, theHitCloner, and theUpdator.

                                                      {
    return new KFTrajectorySmoother(theAlongPropagator,theUpdator,theEstimator,theErrorRescaling,minHits_,theGeometry,theHitCloner);
  }

const MeasurementEstimator* KFTrajectorySmoother::estimator ( void  ) const

inline

Definition at line 84 of file KFTrajectorySmoother.h.

References theEstimator.

Referenced by trajectory().

{return theEstimator;}

const Propagator* KFTrajectorySmoother::oppositePropagator ( ) const

inline

Definition at line 81 of file KFTrajectorySmoother.h.

References theOppositePropagator.

{return theOppositePropagator;}

virtual void KFTrajectorySmoother::setHitCloner ( TkCloner const *  hc )

inlinevirtual

Implements TrajectorySmoother.

Definition at line 91 of file KFTrajectorySmoother.h.

References AnalysisDataFormats_SUSYBSMObjects::hc, and theHitCloner.

{  theHitCloner =   hc;}

Trajectory KFTrajectorySmoother::trajectory ( const Trajectory &  aTraj ) const

overridevirtual

Implements TrajectorySmoother.

Definition at line 36 of file KFTrajectorySmoother.cc.

References alongMomentum, assert(), heavyFlavorValidationHarvestingSequence_cff::combiner, MuonSubdetId::CSC, CSCDetId, DetId::det(), Trajectory::direction(), PXFDetId::disk(), MuonSubdetId::DT, Trajectory::empty(), MeasurementEstimator::estimate(), estimator(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), DetLayerGeometry::idToLayer(), TrajectoryStateOnSurface::isValid(), j, PXBDetId::layer(), TOBDetId::layer(), TIBDetId::layer(), TrajectoryStateOnSurface::localPosition(), LogDebug, LogTrace, TkCloner::makeShared(), Trajectory::measurements(), minHits_, DetId::Muon, nullptr, PV3DBase< T, PVType, FrameType >::perp(), PixelSubdetector::PixelBarrel, PixelSubdetector::PixelEndcap, Propagator::propagate(), Propagator::propagationDirection(), TrajectoryStateOnSurface::rescaleError(), runTheMatrix::ret, MuonSubdetId::RPC, RPCDetId, Trajectory::seed(), dqm_diff::start, DetId::subdetId(), StripSubdetector::TEC, theAlongPropagator, theErrorRescaling, theGeometry, theHitCloner, theOppositePropagator, StripSubdetector::TIB, StripSubdetector::TID, StripSubdetector::TOB, DetId::Tracker, unlikely, TrajectoryStateUpdator::update(), updator(), TIDDetId::wheel(), and TECDetId::wheel().

                                                              {

  if(aTraj.empty()) return Trajectory();

  const Propagator* usePropagator = theAlongPropagator;
  if(aTraj.direction() == alongMomentum) {
    usePropagator = theOppositePropagator;
  }

  const std::vector<TM> & avtm = aTraj.measurements();
  
  
  
  
#ifdef EDM_ML_DEBUG
  LogDebug("TrackFitters") << "KFTrajectorySmoother::trajectories starting with " << avtm.size() << " HITS\n";
  for (unsigned int j=0;j<avtm.size();j++) { 
    if (avtm[j].recHit()->det()) 
      LogTrace("TrackFitters") << "hit #:" << j+1 << " rawId=" << avtm[j].recHit()->det()->geographicalId().rawId() 
                   << " validity=" << avtm[j].recHit()->isValid();
    else
      LogTrace("TrackFitters") << "hit #:" << j+1 << " Hit with no Det information";
  }
#endif // EDM_ML_DEBUG
  
  
  
  TrajectoryStateCombiner combiner;
  bool retry=false;
  auto start = avtm.rbegin();

  do {
    auto hitSize = avtm.rend()-start;
    if unlikely( hitSize < minHits_ ) {
      LogDebug("TrackFitters") << " killing trajectory" << "\n";
      return Trajectory();
    }
    Trajectory ret(aTraj.seed(), usePropagator->propagationDirection());
    Trajectory & myTraj = ret;
    myTraj.reserve(hitSize);
    retry=false;

   TSOS predTsos = (*start).forwardPredictedState();
   predTsos.rescaleError(theErrorRescaling);
   TSOS currTsos;

    auto hitCounter = hitSize;
    for(std::vector<TM>::const_reverse_iterator itm = start; itm != (avtm.rend()); ++itm,--hitCounter) {

    TransientTrackingRecHit::ConstRecHitPointer hit = itm->recHit();

    //check surface just for safety: should never be ==0 because they are skipped in the fitter 
    // if unlikely(hit->det() == nullptr) continue;
    if unlikely( hit->surface()==nullptr ) {
    LogDebug("TrackFitters") << " Error: invalid hit with no GeomDet attached .... skipping";
    continue;
      }


    if (itm != start)//no propagation needed for first smoothed (==last fitted) hit 
      predTsos = usePropagator->propagate( currTsos, *(hit->surface()) );

    if unlikely(!predTsos.isValid()) {
      LogDebug("TrackFitters") << "KFTrajectorySmoother: predicted tsos not valid!";
      LogDebug("TrackFitters") << " retry with last hit removed" << "\n";
      LogDebug("TrackFitters")
      // std::cout 
                << "tsos not valid " << currTsos.globalMomentum().perp() << ' ' 
               << hitSize << ' ' << hitCounter << ' ' << hit->geographicalId().rawId() << ' '
                << hit->surface()->position().perp() << ' ' << hit->surface()->eta() << ' ' << hit->surface()->phi() << std::endl;
      start++;
      retry = true;        
      break;
    }

    if(hit->isValid()) {
 
#ifdef EDM_ML_DEBUG
      LogDebug("TrackFitters")
    << "----------------- HIT #" << hitCounter << " (VALID)-----------------------\n"
    << "HIT IS AT R   " << hit->globalPosition().perp() << "\n"
    << "HIT IS AT Z   " << hit->globalPosition().z() << "\n"
    << "HIT IS AT Phi " << hit->globalPosition().phi() << "\n"
    << "HIT IS AT Loc " << hit->localPosition() << "\n"
    << "WITH LocError " << hit->localPositionError() << "\n"
    << "HIT IS AT Glo " << hit->globalPosition() << "\n"
    << "SURFACE POSITION: " << hit->surface()->position() << "\n"
    << "SURFACE ROTATION: " << hit->surface()->rotation() << "\n"
    << "hit geographicalId=" << hit->geographicalId().rawId();
      
      DetId hitId = hit->geographicalId();
      
      if(hitId.det() == DetId::Tracker) {
    if (hitId.subdetId() == StripSubdetector::TIB )  
      LogTrace("TrackFitters") << " I am TIB " << TIBDetId(hitId).layer();
    else if (hitId.subdetId() == StripSubdetector::TOB ) 
      LogTrace("TrackFitters") << " I am TOB " << TOBDetId(hitId).layer();
    else if (hitId.subdetId() == StripSubdetector::TEC ) 
      LogTrace("TrackFitters") << " I am TEC " << TECDetId(hitId).wheel();
    else if (hitId.subdetId() == StripSubdetector::TID ) 
      LogTrace("TrackFitters") << " I am TID " << TIDDetId(hitId).wheel();
    else if (hitId.subdetId() == (int) PixelSubdetector::PixelBarrel ) 
      LogTrace("TrackFitters") << " I am PixBar " << PXBDetId(hitId).layer();
    else if (hitId.subdetId() == (int) PixelSubdetector::PixelEndcap )
      LogTrace("TrackFitters") << " 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") << " I am DT " << DTWireId(hitId);
    else if (hitId.subdetId() == MuonSubdetId::CSC )
      LogTrace("TrackFitters") << " I am CSC " << CSCDetId(hitId);
    else if (hitId.subdetId() == MuonSubdetId::RPC )
      LogTrace("TrackFitters") << " I am RPC " << RPCDetId(hitId);
    else 
      LogTrace("TrackFitters") << " UNKNOWN MUON HIT TYPE ";
      }
      else
    LogTrace("TrackFitters") << " UNKNOWN HIT TYPE ";
#endif //EDM_ML_DEBUG
      
      
      
      TSOS combTsos,smooTsos;
      
      //3 different possibilities to calculate smoothed state:
      //1: update combined predictions with hit
      //2: combine fwd-prediction with bwd-filter
      //3: combine bwd-prediction with fwd-filter

      //combTsos is the predicted state with N-1 hits information. this means: 
      //forward predicted state for first smoothed (last fitted) hit
      //backward predicted state for last smoothed (first fitted) hit
      //combination of forward and backward predictions for other hits
      if (itm == start) combTsos = itm->forwardPredictedState();
      else if (hitCounter == 1) combTsos = predTsos;
      else combTsos = combiner(predTsos, itm->forwardPredictedState());
      
      if unlikely(!combTsos.isValid()) {
      LogDebug("TrackFitters") << 
        "KFTrajectorySmoother: combined tsos not valid!\n" <<
        "pred Tsos pos: " << predTsos.globalPosition() << "\n" <<
        "pred Tsos mom: " << predTsos.globalMomentum() << "\n" <<
        "TrackingRecHit: " << hit->surface()->toGlobal(hit->localPosition()) << "\n" ;
          start++;
          retry = true;
          break;
    }
      
        assert( (hit->geographicalId()!=0U) | (!hit->canImproveWithTrack()) );
        assert(hit->surface()!=nullptr);
        assert( (!(hit)->canImproveWithTrack()) | (nullptr!=theHitCloner));
        assert( (!(hit)->canImproveWithTrack()) | (nullptr!=dynamic_cast<BaseTrackerRecHit const*>(hit.get())));
        auto preciseHit = theHitCloner->makeShared(hit,combTsos);
        assert(preciseHit->isValid());
        assert( (preciseHit->geographicalId()!=0U) | (!preciseHit->canImproveWithTrack()) );
        assert(preciseHit->surface()!=nullptr);

      if unlikely(!preciseHit->isValid()){
      LogTrace("TrackFitters") << "THE Precise HIT IS NOT VALID: using currTsos = predTsos" << "\n";
      currTsos = predTsos;
      myTraj.push(TM(predTsos, hit, 0, theGeometry->idToLayer(hit->geographicalId()) ));
    }else{
    LogTrace("TrackFitters") << "THE Precise HIT IS VALID: updating currTsos" << "\n";
    
    //update backward predicted tsos with the hit
    currTsos = updator()->update(predTsos, *preciseHit);
        if unlikely(!currTsos.isValid()) {
        currTsos = predTsos;
        edm::LogWarning("KFSmoother_UpdateFailed") << 
          "Failed updating state with hit. Rolling back to non-updated state.\n" <<
          "State: "   << predTsos << 
          "Hit local pos:  " << hit->localPosition() << "\n" <<
          "Hit local err:  " << hit->localPositionError() << "\n" <<
          "Hit global pos: " << hit->globalPosition() << "\n" <<
          "Hit global err: " << hit->globalPositionError().matrix() << 
          "\n";
      }
    
    //smooTsos updates the N-1 hits prediction with the hit
    if (itm == start) smooTsos = itm->updatedState();
    else if (hitCounter == 1) smooTsos = currTsos;
    else smooTsos = combiner(itm->forwardPredictedState(), currTsos); 
    
    if unlikely(!smooTsos.isValid()) {
        LogDebug("TrackFitters") << "KFTrajectorySmoother: smoothed tsos not valid!";
            start++;
            retry = true;
            break;
      }
    
    double estimate;
    if (itm != start) estimate = estimator()->estimate(combTsos, *preciseHit ).second;//correct?
    else estimate = itm->estimate();
    
    LogTrace("TrackFitters")
      << "predTsos !" << "\n"
      << predTsos 
          <<" with local position " << predTsos.localPosition() << "\n\n"
      << "currTsos !" << "\n"
      << currTsos << "\n"
          <<" with local position " << currTsos.localPosition() << "\n\n"
      << "smooTsos !" << "\n"
      << smooTsos 
          <<" with local position " << smooTsos.localPosition() << "\n\n"
      << "smoothing estimate (with combTSOS)=" << estimate << "\n"
      << "filtering estimate=" << itm->estimate() << "\n";
    
    //check for valid hits with no det (refitter with constraints)
    if (preciseHit->det()) myTraj.push(TM(itm->forwardPredictedState(),
                          predTsos,
                          smooTsos,
                          preciseHit,
                          estimate,
                          theGeometry->idToLayer(preciseHit->geographicalId()) ),
                       estimator()->estimate(predTsos,*preciseHit).second);
    else myTraj.push(TM(itm->forwardPredictedState(),
                predTsos,
                smooTsos,
                preciseHit,
                estimate),
             estimator()->estimate(predTsos,*preciseHit).second);
    //itm->estimate());
      }
    } else {
      LogDebug("TrackFitters") 
    << "----------------- HIT #" << hitCounter << " (INVALID)-----------------------";      
      
      //no update
      currTsos = predTsos;
      TSOS combTsos;
      if (itm == start) combTsos = itm->forwardPredictedState();
      else if (hitCounter == 1) combTsos = predTsos;
      else combTsos = combiner(predTsos, itm->forwardPredictedState());
      
      if unlikely(!combTsos.isValid()) {
        LogDebug("TrackFitters") << 
          "KFTrajectorySmoother: combined tsos not valid!";
        return Trajectory();
    }
      assert( (hit->det()==nullptr) || hit->geographicalId()!=0U);
      if (hit->det()) 
         myTraj.push(TM(itm->forwardPredictedState(),
             predTsos,
                 combTsos,
                 hit,
             0,
             theGeometry->idToLayer(hit->geographicalId()) ));
     else myTraj.push(TM(itm->forwardPredictedState(),
                     predTsos,
                     combTsos,
                     hit,
                     0));

    }
  } // for loop

   if (!retry) return ret;
  } while(true);
  
  return Trajectory(); 
  
}

const TrajectoryStateUpdator* KFTrajectorySmoother::updator ( ) const

inline

Definition at line 83 of file KFTrajectorySmoother.h.

References theUpdator.

Referenced by trajectory().

{return theUpdator;}

Member Data Documentation

const DetLayerGeometry KFTrajectorySmoother::dummyGeometry

private

Definition at line 95 of file KFTrajectorySmoother.h.

Referenced by KFTrajectorySmoother().

int KFTrajectorySmoother::minHits_

private

Definition at line 102 of file KFTrajectorySmoother.h.

Referenced by clone(), and trajectory().

const Propagator* KFTrajectorySmoother::theAlongPropagator

private

Definition at line 96 of file KFTrajectorySmoother.h.

Referenced by alongPropagator(), clone(), KFTrajectorySmoother(), trajectory(), and ~KFTrajectorySmoother().

float KFTrajectorySmoother::theErrorRescaling

private

Definition at line 101 of file KFTrajectorySmoother.h.

Referenced by clone(), and trajectory().

const MeasurementEstimator* KFTrajectorySmoother::theEstimator

private

Definition at line 99 of file KFTrajectorySmoother.h.

Referenced by clone(), estimator(), and ~KFTrajectorySmoother().

const DetLayerGeometry* KFTrajectorySmoother::theGeometry

private

Definition at line 103 of file KFTrajectorySmoother.h.

Referenced by clone(), KFTrajectorySmoother(), and trajectory().

TkCloner const* KFTrajectorySmoother::theHitCloner = 0

private

Definition at line 100 of file KFTrajectorySmoother.h.

Referenced by clone(), setHitCloner(), and trajectory().

const Propagator* KFTrajectorySmoother::theOppositePropagator

private

Definition at line 97 of file KFTrajectorySmoother.h.

Referenced by KFTrajectorySmoother(), oppositePropagator(), trajectory(), and ~KFTrajectorySmoother().

const TrajectoryStateUpdator* KFTrajectorySmoother::theUpdator

private

Definition at line 98 of file KFTrajectorySmoother.h.

Referenced by clone(), updator(), and ~KFTrajectorySmoother().