#include <KFFittingSmoother.h>

Inheritance diagram for KFFittingSmoother:

Public Member Functions
std::unique_ptr< TrajectoryFitter >	clone () const override

Trajectory	fitOne (const Trajectory &t, fitType type) const

Trajectory	fitOne (const TrajectorySeed &aSeed, const RecHitContainer &hits, const TrajectoryStateOnSurface &firstPredTsos, fitType type) const

Trajectory	fitOne (const TrajectorySeed &aSeed, const RecHitContainer &hits, fitType type) const

const TrajectoryFitter *	fitter () const

	KFFittingSmoother (const TrajectoryFitter &aFitter, const TrajectorySmoother &aSmoother, double estimateCut=-1, double logPixelProbabilityCut=-16.0, int minNumberOfHits=5, bool rejectTracks=false, bool BreakTrajWith2ConsecutiveMissing=false, bool NoInvalidHitsBeginEnd=false)
	constructor with predefined fitter and smoother and propagator More...

virtual void	setHitCloner (TkCloner const *hc)

const TrajectorySmoother *	smoother () const

virtual	~KFFittingSmoother ()

Public Member Functions inherited from TrajectoryFitter
std::vector< Trajectory >	fit (const Trajectory &traj, fitType type=standard) const

std::vector< Trajectory >	fit (const TrajectorySeed &seed, const RecHitContainer &hits, fitType type=standard) const

std::vector< Trajectory >	fit (const TrajectorySeed &seed, const RecHitContainer &hits, const TrajectoryStateOnSurface &tsos, fitType type=standard) const

virtual	~TrajectoryFitter ()

Private Member Functions
bool	checkForNans (const Trajectory &theTraj) const
	Method to check that the trajectory has no NaN in the states and chi2. More...

Trajectory	smoothingStep (Trajectory const &fitted) const

Private Attributes
bool	breakTrajWith2ConsecutiveMissing

bool	noInvalidHitsBeginEnd

bool	rejectTracksFlag

double	theEstimateCut

const std::unique_ptr < TrajectoryFitter >	theFitter

double	theLogPixelProbabilityCut

int	theMinNumberOfHits

const std::unique_ptr < TrajectorySmoother >	theSmoother

TrajectoryStateWithArbitraryError	tsosWithError

Additional Inherited Members
Public Types inherited from TrajectoryFitter
enum	fitType { standard =0, looper =1 }

typedef Trajectory::RecHitContainer	RecHitContainer

typedef TrajectoryFitterRecord	Record

Detailed Description

A TrajectorySmoother that rpeats the forward fit before smoothing. This is necessary e.g. when the seed introduced a bias (by using a beam contraint etc.). Ported from ORCA

Author: todorov, cerati

Definition at line 17 of file KFFittingSmoother.h.

Constructor & Destructor Documentation

KFFittingSmoother::KFFittingSmoother	(	const TrajectoryFitter &	aFitter,
		const TrajectorySmoother &	aSmoother,
		double	estimateCut = `-1`,
		double	logPixelProbabilityCut = `-16.0`,
		int	minNumberOfHits = `5`,
		bool	rejectTracks = `false`,
		bool	BreakTrajWith2ConsecutiveMissing = `false`,
		bool	NoInvalidHitsBeginEnd = `false`
	)

inline

constructor with predefined fitter and smoother and propagator

Definition at line 21 of file KFFittingSmoother.h.

Referenced by clone().

                                                       :
   theFitter(aFitter.clone()),
       theSmoother(aSmoother.clone()),
       theEstimateCut(estimateCut),
 
       // ggiurgiu@fnal.gov
       theLogPixelProbabilityCut( logPixelProbabilityCut ),
 
       theMinNumberOfHits(minNumberOfHits),
       rejectTracksFlag(rejectTracks),
       breakTrajWith2ConsecutiveMissing(BreakTrajWith2ConsecutiveMissing),
       noInvalidHitsBeginEnd(NoInvalidHitsBeginEnd) {}

virtual KFFittingSmoother::~KFFittingSmoother ( )

inlinevirtual

Definition at line 41 of file KFFittingSmoother.h.

41 {}

Member Function Documentation

bool KFFittingSmoother::checkForNans ( const Trajectory & theTraj ) const

private

Method to check that the trajectory has no NaN in the states and chi2.

Definition at line 22 of file KFFittingSmoother.cc.

References Trajectory::chiSquared(), i, edm::isNotFinite(), j, visualization-live-secondInstance_cfg::m, Trajectory::measurements(), and findQualityFiles::v.

 {
   if (edm::isNotFinite(theTraj.chiSquared ())) return false;
   const std::vector<TrajectoryMeasurement> & vtm = theTraj.measurements();
   for (std::vector<TrajectoryMeasurement>::const_iterator tm=vtm.begin(); tm!=vtm.end();tm++) {
     if (edm::isNotFinite(tm->estimate())) return false;
     AlgebraicVector5 v = tm->updatedState ().localParameters ().vector();
     for (int i=0;i<5;++i) if (edm::isNotFinite(v[i])) return false;
     const AlgebraicSymMatrix55 & m = tm->updatedState ().curvilinearError ().matrix();
     for (int i=0;i<5;++i)
       for (int j=0;j<(i+1);++j) if (edm::isNotFinite(m(i,j))) return false;
   }
   return true;
 }

std::unique_ptr<TrajectoryFitter> KFFittingSmoother::clone ( void ) const

inlineoverridevirtual

Implements TrajectoryFitter.

Definition at line 53 of file KFFittingSmoother.h.

References breakTrajWith2ConsecutiveMissing, KFFittingSmoother(), noInvalidHitsBeginEnd, rejectTracksFlag, theEstimateCut, theFitter, theLogPixelProbabilityCut, theMinNumberOfHits, and theSmoother.

                                                          {
     return std::unique_ptr<TrajectoryFitter>(
         new KFFittingSmoother(*theFitter,
                               *theSmoother,
                               theEstimateCut,theLogPixelProbabilityCut,
                               theMinNumberOfHits,rejectTracksFlag,
                               breakTrajWith2ConsecutiveMissing,noInvalidHitsBeginEnd));
   }

Trajectory KFFittingSmoother::fitOne	(	const Trajectory &	t,
		fitType	type
	)		const

virtual

Implements TrajectoryFitter.

Definition at line 17 of file KFFittingSmoother.cc.

References Trajectory::isValid().

                                                                             {
   if (!t.isValid() ) return Trajectory();
   return smoothingStep(theFitter->fitOne(t,type));
 }

Trajectory KFFittingSmoother::fitOne	(	const TrajectorySeed &	aSeed,
		const RecHitContainer &	hits,
		const TrajectoryStateOnSurface &	firstPredTsos,
		fitType	type
	)		const

virtual

Implements TrajectoryFitter.

Definition at line 37 of file KFFittingSmoother.cc.

References Trajectory::chiSquared(), SiPixelRecHit::clusterProbability(), GOODCOLL_filter_cfg::cut, Trajectory::direction(), Trajectory::firstMeasurement(), Trajectory::foundHits(), TrackingRecHit::hit(), TrackingRecHit::isValid(), Trajectory::isValid(), j, Trajectory::lastMeasurement(), LogDebug, LogTrace, Trajectory::measurements(), TrackingRecHit::missing, eostools::move(), PixelSubdetector::PixelBarrel, PixelSubdetector::PixelEndcap, Trajectory::pop(), SiPixelRecHit::qBin(), TrajectoryMeasurement::recHitR(), Trajectory::seed(), std::swap(), and DetId::Tracker.

 {
   LogDebug("TrackFitters") << "In KFFittingSmoother::fit";
 
   if ( hits.empty() ) return Trajectory();
 
   bool hasoutliers=false;
   bool has_low_pixel_prob=false; // ggiurgiu@fnal.gov: Add flag for pixel hits with low template probability
 
   // ggiurgiu@fnal.gov: If log(Prob) < -15.0 or if Prob <= 0.0 then set log(Prob) = -15.0
   double log_pixel_probability_lower_limit = -15.0;
 
   RecHitContainer myHits = hits;
   Trajectory smoothed;
   Trajectory tmp_first;
   bool firstTry=true;
 
   do {
 
     //if no outliers the fit is done only once
     //for (unsigned int j=0;j<myHits.size();j++) {
     //if (myHits[j]->det())
     //LogTrace("TrackFitters") << "hit #:" << j+1 << " rawId=" << myHits[j]->det()->geographicalId().rawId()
     //<< " validity=" << myHits[j]->isValid();
     //else
     //LogTrace("TrackFitters") << "hit #:" << j+1 << " Hit with no Det information";
     //}
 
     hasoutliers        = false;
     has_low_pixel_prob = false; // ggiurgiu@fnal.gov
 
     double cut = theEstimateCut;
 
     double log_pixel_prob_cut = theLogPixelProbabilityCut;  // ggiurgiu@fnal.gov
 
 
     unsigned int outlierId = 0xffffffff;
     const GeomDet* outlierDet = 0;
 
     unsigned int low_pixel_prob_Id = 0; // ggiurgiu@fnal.gov
     const GeomDet* low_pixel_prob_Det = 0; // ggiurgiu@fnal.gov
 
     //call the fitter
     smoothed  = smoothingStep(theFitter->fitOne(aSeed, myHits, firstPredTsos));
 
     //if (tmp_first.size()==0) tmp_first = smoothed; moved later
 
     if ( !smoothed.isValid() )  {
       if ( rejectTracksFlag ) {
     LogTrace("TrackFitters") << "smoothed invalid => trajectory rejected";
     //return vector<Trajectory>(); // break is enough to get this
       } else {
     LogTrace("TrackFitters") << "smoothed invalid => returning orignal trajectory" ;
     std::swap(smoothed, tmp_first); // if first attempt, tmp_first would be invalid anyway
       }
       break;
     }
     //else {
     //LogTrace("TrackFitters") << "dump hits after smoothing";
     //Trajectory::DataContainer meas = smoothed[0].measurements();
     //for (Trajectory::DataContainer::iterator it=meas.begin();it!=meas.end();++it) {
     //LogTrace("TrackFitters") << "hit #" << meas.end()-it-1 << " validity=" << it->recHit()->isValid()
     //<< " det=" << it->recHit()->geographicalId().rawId();
     //}
     //}
 
     if (!checkForNans(smoothed)) {
       edm::LogWarning("TrackNaN")<<"Track has NaN";
       return Trajectory();
     }
 
 
     if ( theEstimateCut > 0 || log_pixel_prob_cut > log_pixel_probability_lower_limit ) {
       if ( smoothed.foundHits() < theMinNumberOfHits ) {
     if ( rejectTracksFlag ) {
       LogTrace("TrackFitters") << "smoothed.foundHits()<theMinNumberOfHits => trajectory rejected";
       return Trajectory(); // invalid
     } else{
       // it might be it's the first step
       if ( firstTry ) {
         firstTry = false;
         std::swap(tmp_first,smoothed);
       }
 
       LogTrace("TrackFitters")
         << "smoothed.foundHits()<theMinNumberOfHits => returning orignal trajectory with chi2="
         <<  smoothed.chiSquared() ;
     }
     break;
       }
 
       // Check if there are outliers or low probability pixel rec hits
       const std::vector<TrajectoryMeasurement> & vtm = smoothed.measurements();
 
       double log_pixel_hit_probability = -999999.9;
 
       for (std::vector<TrajectoryMeasurement>::const_iterator tm=vtm.begin(); tm!=vtm.end();tm++) {
     double estimate = tm->estimate();
 
     // --- here is the block of code about generic chi2-based Outlier Rejection ---
     if ( ( (estimate>cut) & (theEstimateCut>0) ) && 
              tm->recHit()->det()!=nullptr // do not consider outliers constraints and other special "hits" 
            ) {  
       hasoutliers = true;
       cut = estimate;
       outlierId  = tm->recHit()->geographicalId().rawId();
       outlierDet = tm->recHit()->det();
     }
     // --- here the block of code about generic chi2-based Outlier Rejection ends ---
 
 
     // --- here is the block of code about PXL Outlier Rejection ---
     if(log_pixel_prob_cut > log_pixel_probability_lower_limit){
       // TO BE FIXED: the following code should really be moved into an external class or
       // at least in a separate function. Current solution is ugly!
       // The KFFittingSmoother shouldn't handle the details of
       // Outliers identification and rejection. It shoudl just fit tracks.
 
       // ggiurgiu@fnal.gov: Get pixel hit probability here
       TransientTrackingRecHit::ConstRecHitPointer hit = tm->recHit();
       unsigned int testSubDetID = ( hit->geographicalId().subdetId() );
 
       if ( hit->isValid() &&
            hit->geographicalId().det() == DetId::Tracker &&
            ( testSubDetID == PixelSubdetector::PixelBarrel ||
          testSubDetID == PixelSubdetector::PixelEndcap )
            ){
         // get the enclosed persistent hit
         const TrackingRecHit* persistentHit = hit->hit();
 
         // check if it's not null, and if it's a valid pixel hit
         if ( !persistentHit == 0 &&
          typeid(*persistentHit) == typeid(SiPixelRecHit) )
           {
 
         // tell the C++ compiler that the hit is a pixel hit
         const SiPixelRecHit* pixhit = dynamic_cast<const SiPixelRecHit*>( hit->hit() );
 
         double pixel_hit_probability = (float)pixhit->clusterProbability(0);
 
         if ( pixel_hit_probability < 0.0 )
           LogDebug("From KFFittingSmoother.cc")
             << "Wraning : Negative pixel hit probability !!!! Will set the probability to 10^{-15} !!!" << endl;
 
         if ( pixel_hit_probability <= 0.0 || log10( pixel_hit_probability ) < log_pixel_probability_lower_limit )
           log_pixel_hit_probability = log_pixel_probability_lower_limit;
         else
           log_pixel_hit_probability = log10( pixel_hit_probability );
 
         int qbin = (int)pixhit->qBin();
 
         if ( ( log_pixel_hit_probability <  log_pixel_prob_cut ) &&
              ( qbin                      != 0                  ) )
           {
             has_low_pixel_prob = true;
             log_pixel_prob_cut = log_pixel_hit_probability;
             low_pixel_prob_Id  = tm->recHit()->geographicalId().rawId();
             low_pixel_prob_Det = tm->recHit()->det();
           }
 
           } // if ( !persistentHit == 0 && ... )
 
       } // if ( hit->isValid() && ... )
     }
     // --- here the block of code about PXL Outlier Rejection ends ---
 
       } // for (std::vector<TrajectoryMeasurement>::const_iterator tm=vtm.begin(); tm!=vtm.end(); tm++)
 
 
       if ( hasoutliers || has_low_pixel_prob ) { // Reject outliers and pixel hits with low probability
 
     // Replace outlier hit or low probability pixel hit with invalid hit
     for ( unsigned int j=0; j<myHits.size(); ++j ) {
       if ( hasoutliers && outlierId == myHits[j]->geographicalId().rawId() )
         {
           LogTrace("TrackFitters") << "Rejecting outlier hit  with estimate " << cut << " at position "
                        << j << " with rawId=" << myHits[j]->geographicalId().rawId();
           LogTrace("TrackFitters") << "The fit will be repeated without the outlier";
           myHits[j] = std::make_shared<InvalidTrackingRecHit>(*outlierDet, TrackingRecHit::missing);
         }
       else if ( has_low_pixel_prob && low_pixel_prob_Id == myHits[j]->geographicalId().rawId() ){
         LogTrace("TrackFitters") << "Rejecting low proability pixel hit with log_pixel_prob_cut = "
                      << log_pixel_prob_cut << " at position "
                      << j << " with rawId =" << myHits[j]->geographicalId().rawId();
         LogTrace("TrackFitters") << "The fit will be repeated without the outlier";
         myHits[j] = std::make_shared<InvalidTrackingRecHit>(*low_pixel_prob_Det, TrackingRecHit::missing);
       }
 
     } // for ( unsigned int j=0; j<myHits.size(); ++j)
 
     // Look if there are two consecutive invalid hits
     if ( breakTrajWith2ConsecutiveMissing ) {
       unsigned int firstinvalid = myHits.size();
       for ( unsigned int j=0; j<myHits.size()-1; ++j )
         {
           if ( ((myHits[j  ]->type() == TrackingRecHit::missing) && (myHits[j  ]->geographicalId().rawId() != 0)) &&
            ((myHits[j+1]->type() == TrackingRecHit::missing) && (myHits[j+1]->geographicalId().rawId() != 0)) )
         {
           firstinvalid = j;
           LogTrace("TrackFitters") << "Found two consecutive missing hits. First invalid: " << firstinvalid;
           break;
         }
         }
 
       //reject all the hits after the last valid before two consecutive invalid (missing) hits
       //hits are sorted in the same order as in the track candidate FIXME??????
       if (firstinvalid != myHits.size()) myHits.erase(myHits.begin()+firstinvalid,myHits.end());
 
     }
 
       } // if ( hasoutliers || has_low_pixel_prob )
 
     } // if ( theEstimateCut > 0 ... )
 
     if ( ( hasoutliers ||        // otherwise there won't be a 'next' loop where tmp_first could be useful
        has_low_pixel_prob ) &&  // ggiurgiu@fnal.gov
      !rejectTracksFlag &&  // othewrise we won't ever need tmp_first
      firstTry ) { // only at first step
       std::swap(smoothed,tmp_first); firstTry=false;
     }
 
   } // do
   while ( hasoutliers || has_low_pixel_prob ); // ggiurgiu@fnal.gov
 
   if ( smoothed.isValid() ) {
     if ( noInvalidHitsBeginEnd )  {
       // discard latest dummy measurements
       if (!smoothed.lastMeasurement().recHitR().isValid() )
     LogTrace("TrackFitters") << "Last measurement is invalid";
 
       while (!smoothed.lastMeasurement().recHitR().isValid() )
     smoothed.pop();
 
       //remove the invalid hits at the begin of the trajectory
       if ( !smoothed.firstMeasurement().recHitR().isValid() ) {
     LogTrace("TrackFitters") << "First measurement is in`valid";
     Trajectory tmpTraj(smoothed.seed(),smoothed.direction());
     Trajectory::DataContainer  & meas = smoothed.measurements();
     auto it = meas.begin();
     for ( ; it!=meas.end(); ++it )
       if ( it->recHitR().isValid() )  break;
     tmpTraj.push(std::move(*it),smoothed.chiSquared());//push the first valid measurement and set the same global chi2
 
     for (auto itt=it+1; itt!=meas.end();++itt)
       tmpTraj.push(std::move(*itt),0);//add all the other measurements
 
     std::swap(smoothed,tmpTraj);
 
       } //  if ( !smoothed[0].firstMeasurement().recHit()->isValid() )
 
     } // if ( noInvalidHitsBeginEnd )
 
     LogTrace("TrackFitters") << "end: returning smoothed trajectory with chi2="
                  << smoothed.chiSquared() ;
 
     //LogTrace("TrackFitters") << "dump hits before return";
     //Trajectory::DataContainer meas = smoothed[0].measurements();
     //for (Trajectory::DataContainer::iterator it=meas.begin();it!=meas.end();++it) {
     //LogTrace("TrackFitters") << "hit #" << meas.end()-it-1 << " validity=" << it->recHit()->isValid()
     //<< " det=" << it->recHit()->geographicalId().rawId();
     //}
 
   }
 
   return smoothed;
 
 }

Trajectory KFFittingSmoother::fitOne	(	const TrajectorySeed &	aSeed,
		const RecHitContainer &	hits,
		fitType	type
	)		const

virtual

Implements TrajectoryFitter.

Definition at line 309 of file KFFittingSmoother.cc.

References edm::hlt::Exception.

                                                                     {
 
   throw cms::Exception("TrackFitters",
                "KFFittingSmoother::fit(TrajectorySeed, <TransientTrackingRecHit>) not implemented");
 
   return Trajectory();
 }

const TrajectoryFitter* KFFittingSmoother::fitter ( void ) const

inline

Definition at line 50 of file KFFittingSmoother.h.

References theFitter.

Referenced by KalmanAlignmentAlgorithm::initializeAlignmentSetups().

50 {return theFitter.get();}

KFFittingSmoother::theFitter

const std::unique_ptr< TrajectoryFitter > theFitter

Definition: KFFittingSmoother.h:75

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

inlinevirtual

Implements TrajectoryFitter.

Definition at line 63 of file KFFittingSmoother.h.

References theFitter, and theSmoother.

                                                  {
         theFitter->setHitCloner(hc);
         theSmoother->setHitCloner(hc);
   }

const TrajectorySmoother* KFFittingSmoother::smoother ( ) const

inline

Definition at line 51 of file KFFittingSmoother.h.

References theSmoother.

Referenced by KalmanAlignmentAlgorithm::initializeAlignmentSetups().

51 {return theSmoother.get();}

KFFittingSmoother::theSmoother

const std::unique_ptr< TrajectorySmoother > theSmoother

Definition: KFFittingSmoother.h:76

Trajectory KFFittingSmoother::smoothingStep ( Trajectory const & fitted ) const

inlineprivate

Definition at line 71 of file KFFittingSmoother.h.

References theSmoother.

71 {return theSmoother->trajectory(fitted);}

KFFittingSmoother::theSmoother

const std::unique_ptr< TrajectorySmoother > theSmoother

Definition: KFFittingSmoother.h:76

Member Data Documentation

bool KFFittingSmoother::breakTrajWith2ConsecutiveMissing

private

Definition at line 83 of file KFFittingSmoother.h.

Referenced by clone().

bool KFFittingSmoother::noInvalidHitsBeginEnd

private

Definition at line 84 of file KFFittingSmoother.h.

Referenced by clone().

bool KFFittingSmoother::rejectTracksFlag

private

Definition at line 82 of file KFFittingSmoother.h.

Referenced by clone().

double KFFittingSmoother::theEstimateCut

private

Definition at line 77 of file KFFittingSmoother.h.

Referenced by clone().

const std::unique_ptr<TrajectoryFitter> KFFittingSmoother::theFitter

private

Definition at line 75 of file KFFittingSmoother.h.

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

double KFFittingSmoother::theLogPixelProbabilityCut

private

Definition at line 79 of file KFFittingSmoother.h.

Referenced by clone().

int KFFittingSmoother::theMinNumberOfHits

private

Definition at line 81 of file KFFittingSmoother.h.

Referenced by clone().

const std::unique_ptr<TrajectorySmoother> KFFittingSmoother::theSmoother

private

Definition at line 76 of file KFFittingSmoother.h.

Referenced by clone(), setHitCloner(), smoother(), and smoothingStep().

TrajectoryStateWithArbitraryError KFFittingSmoother::tsosWithError

private

Definition at line 86 of file KFFittingSmoother.h.

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation