#include <KFTrajectorySmoother.h>

Inheritance diagram for KFTrajectorySmoother:

Public Member Functions
const Propagator *	alongPropagator () const

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

void	setHitCloner (TkCloner const *hc) override

Trajectory	trajectory (const Trajectory &aTraj) const override

const TrajectoryStateUpdator *	updator () const

	~KFTrajectorySmoother () override

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, theOppositePropagator, trajectory(), and ~KFTrajectorySmoother().

                                                     :
     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 ( )

override

Definition at line 14 of file KFTrajectorySmoother.cc.

References theAlongPropagator, theEstimator, theOppositePropagator, and theUpdator.

Referenced by KFTrajectorySmoother().

                                             {
 
   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.

80 { return theAlongPropagator;}

KFTrajectorySmoother::theAlongPropagator

const Propagator * theAlongPropagator

Definition: KFTrajectorySmoother.h:96

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().

84 {return theEstimator;}

KFTrajectorySmoother::theEstimator

const MeasurementEstimator * theEstimator

Definition: KFTrajectorySmoother.h:99

const Propagator* KFTrajectorySmoother::oppositePropagator ( ) const

inline

Definition at line 81 of file KFTrajectorySmoother.h.

References theOppositePropagator.

81 {return theOppositePropagator;}

KFTrajectorySmoother::theOppositePropagator

const Propagator * theOppositePropagator

Definition: KFTrajectorySmoother.h:97

void KFTrajectorySmoother::setHitCloner ( TkCloner const * hc )

inlineoverridevirtual

Implements TrajectorySmoother.

Definition at line 91 of file KFTrajectorySmoother.h.

References AnalysisDataFormats_SUSYBSMObjects::hc, and theHitCloner.

91 { theHitCloner = hc;}

AnalysisDataFormats_SUSYBSMObjects::hc

susybsm::HSCParticleCollection hc

Definition: classes.h:25

KFTrajectorySmoother::theHitCloner

TkCloner const * theHitCloner

Definition: KFTrajectorySmoother.h:100

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

overridevirtual

Implements TrajectorySmoother.

Definition at line 24 of file KFTrajectorySmoother.cc.

References alongMomentum, Trajectory::direction(), FrontierConditions_GlobalTag_cff::dump, Trajectory::empty(), MeasurementEstimator::estimate(), estimator(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), DetLayerGeometry::idToLayer(), createfilelist::int, TrajectoryStateOnSurface::isValid(), TrajectoryStateOnSurface::localPosition(), LogDebug, LogTrace, TkCloner::makeShared(), Trajectory::measurements(), minHits_, nullptr, PV3DBase< T, PVType, FrameType >::perp(), Propagator::propagate(), Propagator::propagationDirection(), rpcPointValidation_cfi::recHit, TrajectoryStateOnSurface::rescaleError(), Trajectory::reserve(), Trajectory::seed(), command_line::start, theAlongPropagator, theErrorRescaling, theGeometry, theHitCloner, theOppositePropagator, mitigatedMETSequence_cff::U, UNLIKELY, TrajectoryStateUpdator::update(), and updator().

Referenced by KFTrajectorySmoother().

                                                               {
 
   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 << ' ' << int(hit->geographicalId()) << ' '  
                 << hit->surface()->position().perp() << ' ' << hit->surface()->eta() << ' ' << hit->surface()->phi() << std::endl;
       start++;
       retry = true;        
       break;
     }
 
     if(hit->isValid()) {
  
       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);
 
         dump(*hit,hitCounter,"TrackFitters");
 
       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(); 
   
 }