#include <KFFittingSmoother.h>

Inheritance diagram for KFFittingSmoother:

Public Member Functions
KFFittingSmoother *	clone () const

virtual std::vector< Trajectory >	fit (const Trajectory &t) const

virtual std::vector< Trajectory >	fit (const TrajectorySeed &aSeed, const RecHitContainer &hits, const TrajectoryStateOnSurface &firstPredTsos) const

virtual std::vector< Trajectory >	fit (const TrajectorySeed &aSeed, const RecHitContainer &hits) 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...

const TrajectorySmoother *	smoother () const

virtual	~KFFittingSmoother ()

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

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

virtual std::vector< Trajectory >	fit (const TrajectorySeed &seed, const RecHitContainer &hits, const TrajectoryStateOnSurface &tsos, fitType type) 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...

void	smoothingStep (std::vector< Trajectory > &fitted, std::vector< Trajectory > &smoothed) const

Private Attributes
bool	breakTrajWith2ConsecutiveMissing

bool	noInvalidHitsBeginEnd

bool	rejectTracksFlag

double	theEstimateCut

const TrajectoryFitter *	theFitter

double	theLogPixelProbabilityCut

int	theMinNumberOfHits

const 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

Date:: 2012/09/13 08:48:07

Revision:: 1.16.2.1

Author: todorov, cerati

Definition at line 19 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 23 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) {}

KFFittingSmoother::~KFFittingSmoother ( )

virtual

Definition at line 18 of file KFFittingSmoother.cc.

 {
   delete theSmoother;
   delete theFitter;
 }

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 35 of file KFFittingSmoother.cc.

References Trajectory::chiSquared(), i, edm::detail::isnan(), j, m, Trajectory::measurements(), and v.

 {
   if (std::isnan(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 (std::isnan(tm->estimate())) return false;
     AlgebraicVector5 v = tm->updatedState ().localParameters ().vector();
     for (int i=0;i<5;++i) if (std::isnan(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 (std::isnan(m(i,j))) return false;
   }
   return true;
 }

KFFittingSmoother* KFFittingSmoother::clone ( void ) const

inlinevirtual

Implements TrajectoryFitter.

Definition at line 55 of file KFFittingSmoother.h.

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

                                    {
     return new KFFittingSmoother(*theFitter,*theSmoother,
                  theEstimateCut,theLogPixelProbabilityCut,
                  theMinNumberOfHits,rejectTracksFlag,
                  breakTrajWith2ConsecutiveMissing,noInvalidHitsBeginEnd);
   }

vector< Trajectory > KFFittingSmoother::fit ( const Trajectory & t ) const

virtual

Implements TrajectoryFitter.

Definition at line 24 of file KFFittingSmoother.cc.

References Trajectory::isValid().

 {
   vector<Trajectory> smoothed;
   if ( t.isValid() ) 
     { 
       vector<Trajectory> fitted = fitter()->fit(t);
       smoothingStep(fitted, smoothed);
     }
   return smoothed;
 }

vector< Trajectory > KFFittingSmoother::fit	(	const TrajectorySeed &	aSeed,
		const RecHitContainer &	hits,
		const TrajectoryStateOnSurface &	firstPredTsos
	)		const

virtual

Implements TrajectoryFitter.

Definition at line 50 of file KFFittingSmoother.cc.

References newFWLiteAna::build, GOODCOLL_filter_cfg::cut, relativeConstraints::empty, j, LogDebug, LogTrace, TrackingRecHit::missing, PixelSubdetector::PixelBarrel, PixelSubdetector::PixelEndcap, python.multivaluedict::pop(), Trajectory::push(), and DetId::Tracker.

 {
   LogDebug("TrackFitters") << "In KFFittingSmoother::fit";
 
   //if(hits.empty()) return vector<Trajectory>(); // gio: moved later to optimize return value
   
   bool hasoutliers;
   bool has_low_pixel_prob; // 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; 
   vector<Trajectory> smoothed;
   vector<Trajectory> tmp_first;
 
   do
     {
       if ( hits.empty() ) 
     { 
       smoothed.clear(); 
       break; 
     }
       
       //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 = 0;
       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
       vector<Trajectory> fitted = fitter()->fit(aSeed, myHits, firstPredTsos);
       //call the smoother
       smoothed.clear();
       smoothingStep(fitted, smoothed);
       
       //if (tmp_first.size()==0) tmp_first = smoothed; moved later
       
       if ( smoothed.empty() ) 
     {
       if ( rejectTracksFlag )
         {
           LogTrace("TrackFitters") << "smoothed.size()==0 => trajectory rejected";
           //return vector<Trajectory>(); // break is enough to get this
         } 
       else 
         {
           LogTrace("TrackFitters") << "smoothed.size()==0 => returning orignal trajectory" ;
           smoothed.swap(tmp_first); // if first attempt, tmp_first would be empty 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[0])) {
     edm::LogError("TrackNaN")<<"Track has NaN";
     smoothed.clear();
     break;
       }
       
       if ( theEstimateCut > 0 || log_pixel_prob_cut > log_pixel_probability_lower_limit ) 
     {
       if ( smoothed[0].foundHits() < theMinNumberOfHits ) 
         {
           if ( rejectTracksFlag ) 
         {
           LogTrace("TrackFitters") << "smoothed[0].foundHits()<theMinNumberOfHits => trajectory rejected";
           smoothed.clear();
           //return vector<Trajectory>(); // break is enough
         } 
           else 
         {
           // it might be it's the first step
           if ( !tmp_first.empty() ) 
             { 
               tmp_first.swap(smoothed); 
             } 
           
           LogTrace("TrackFitters") 
             << "smoothed[0].foundHits()<theMinNumberOfHits => returning orignal trajectory with chi2=" 
             <<  smoothed[0].chiSquared() ;
         }
           break;
         }
       
       // Check if there are outliers or low probability pixel rec hits
       const std::vector<TrajectoryMeasurement> & vtm = smoothed[0].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 ) 
         {
           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] = InvalidTransientRecHit::build(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] = InvalidTransientRecHit::build(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()-1;
           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??????
           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
        tmp_first.empty() ) 
     { // only at first step
       smoothed.swap(tmp_first);
     }   
       
     } // do
   while ( hasoutliers || has_low_pixel_prob ); // ggiurgiu@fnal.gov
   
   if ( !smoothed.empty() ) 
     {
       if ( noInvalidHitsBeginEnd ) 
     {
       // discard latest dummy measurements
       if ( !smoothed[0].empty() && 
            !smoothed[0].lastMeasurement().recHit()->isValid() ) 
         LogTrace("TrackFitters") << "Last measurement is invalid";
     
       while ( !smoothed[0].empty() && 
           !smoothed[0].lastMeasurement().recHit()->isValid() ) 
         smoothed[0].pop();
       
       //remove the invalid hits at the begin of the trajectory
       if ( !smoothed[0].firstMeasurement().recHit()->isValid() ) 
         {
           LogTrace("TrackFitters") << "First measurement is invalid";
           Trajectory tmpTraj(smoothed[0].seed(),smoothed[0].direction());
           Trajectory::DataContainer meas = smoothed[0].measurements();
           
           Trajectory::DataContainer::iterator it;//first valid hit
           for ( it=meas.begin(); it!=meas.end(); ++it ) 
         {
           if ( !it->recHit()->isValid() ) 
             continue;
           else break;
         }
           tmpTraj.push(*it,smoothed[0].chiSquared());//push the first valid measurement and set the same global chi2
          
           for (Trajectory::DataContainer::iterator itt=it+1; itt!=meas.end();++itt) 
         {
           tmpTraj.push(*itt,0);//add all the other measurements
         }
           
           smoothed.clear();
           smoothed.push_back(tmpTraj);
        
         } //  if ( !smoothed[0].firstMeasurement().recHit()->isValid() ) 
     
     } // if ( noInvalidHitsBeginEnd ) 
       
       LogTrace("TrackFitters") << "end: returning smoothed trajectory with chi2=" 
                    << smoothed[0].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;
 
 }

vector< Trajectory > KFFittingSmoother::fit	(	const TrajectorySeed &	aSeed,
		const RecHitContainer &	hits
	)		const

virtual

Implements TrajectoryFitter.

Definition at line 370 of file KFFittingSmoother.cc.

References edm::hlt::Exception.

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

const TrajectoryFitter* KFFittingSmoother::fitter ( void ) const

inline

Definition at line 52 of file KFFittingSmoother.h.

References theFitter.

Referenced by KalmanAlignmentAlgorithm::initializeAlignmentSetups().

52 {return theFitter;}

KFFittingSmoother::theFitter

const TrajectoryFitter * theFitter

Definition: KFFittingSmoother.h:64

const TrajectorySmoother* KFFittingSmoother::smoother ( ) const

inline

Definition at line 53 of file KFFittingSmoother.h.

References theSmoother.

Referenced by KalmanAlignmentAlgorithm::initializeAlignmentSetups().

53 {return theSmoother;}

KFFittingSmoother::theSmoother

const TrajectorySmoother * theSmoother

Definition: KFFittingSmoother.h:65

void KFFittingSmoother::smoothingStep	(	std::vector< Trajectory > &	fitted,
		std::vector< Trajectory > &	smoothed
	)		const

private

Definition at line 359 of file KFFittingSmoother.cc.

References LogDebug.

 {
  
   for(vector<Trajectory>::iterator it = fitted.begin(); it != fitted.end();
       it++) {
     vector<Trajectory> mysmoothed = smoother()->trajectories(*it);
     smoothed.insert(smoothed.end(), mysmoothed.begin(), mysmoothed.end());
   }
   LogDebug("TrackFitters") << "In KFFittingSmoother::smoothingStep "<<smoothed.size();
 }

Member Data Documentation

bool KFFittingSmoother::breakTrajWith2ConsecutiveMissing

private

Definition at line 72 of file KFFittingSmoother.h.

Referenced by clone().

bool KFFittingSmoother::noInvalidHitsBeginEnd

private

Definition at line 73 of file KFFittingSmoother.h.

Referenced by clone().

bool KFFittingSmoother::rejectTracksFlag

private

Definition at line 71 of file KFFittingSmoother.h.

Referenced by clone().

double KFFittingSmoother::theEstimateCut

private

Definition at line 66 of file KFFittingSmoother.h.

Referenced by clone().

const TrajectoryFitter* KFFittingSmoother::theFitter

private

Definition at line 64 of file KFFittingSmoother.h.

Referenced by clone(), and fitter().

double KFFittingSmoother::theLogPixelProbabilityCut

private

Definition at line 68 of file KFFittingSmoother.h.

Referenced by clone().

int KFFittingSmoother::theMinNumberOfHits

private

Definition at line 70 of file KFFittingSmoother.h.

Referenced by clone().

const TrajectorySmoother* KFFittingSmoother::theSmoother

private

Definition at line 65 of file KFFittingSmoother.h.

Referenced by clone(), and smoother().

TrajectoryStateWithArbitraryError KFFittingSmoother::tsosWithError

private

Definition at line 76 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