CMS 3D CMS Logo

Public Member Functions | Private Member Functions | Private Attributes

KFFittingSmoother Class Reference

#include <KFFittingSmoother.h>

Inheritance diagram for KFFittingSmoother:
TrajectoryFitter

List of all members.

Public Member Functions

KFFittingSmootherclone () const
virtual std::vector< Trajectoryfit (const Trajectory &t) const
virtual std::vector< Trajectoryfit (const TrajectorySeed &aSeed, const RecHitContainer &hits, const TrajectoryStateOnSurface &firstPredTsos) const
virtual std::vector< Trajectoryfit (const TrajectorySeed &aSeed, const RecHitContainer &hits) const
const TrajectoryFitterfitter () 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
const TrajectorySmoothersmoother () const
virtual ~KFFittingSmoother ()

Private Member Functions

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

Private Attributes

bool breakTrajWith2ConsecutiveMissing
bool noInvalidHitsBeginEnd
bool rejectTracksFlag
double theEstimateCut
const TrajectoryFittertheFitter
double theLogPixelProbabilityCut
int theMinNumberOfHits
const TrajectorySmoothertheSmoother
TrajectoryStateWithArbitraryError tsosWithError

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:
2011/01/06 17:04:16
Revision:
1.16
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 17 of file KFFittingSmoother.cc.

{
  delete theSmoother;
  delete theFitter;
}

Member Function Documentation

KFFittingSmoother* KFFittingSmoother::clone ( void  ) const [inline, virtual]
vector< Trajectory > KFFittingSmoother::fit ( const TrajectorySeed aSeed,
const RecHitContainer hits,
const TrajectoryStateOnSurface firstPredTsos 
) const [virtual]

Implements TrajectoryFitter.

Definition at line 34 of file KFFittingSmoother.cc.

References newFWLiteAna::build, SiPixelRecHit::clusterProbability(), align_tpl::cut, relativeConstraints::empty, j, LogDebug, LogTrace, TrackingRecHit::missing, combine::missing, PixelSubdetector::PixelBarrel, GeomDetEnumerators::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 ( 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 ) 
                {
                  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 );
                        
                        if ( log_pixel_hit_probability < log_pixel_prob_cut )
                          {
                            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 346 of file KFFittingSmoother.cc.

References Exception.

                                                                            {

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

  return vector<Trajectory>();
}
vector< Trajectory > KFFittingSmoother::fit ( const Trajectory t) const [virtual]

Implements TrajectoryFitter.

Definition at line 23 of file KFFittingSmoother.cc.

References Trajectory::isValid().

{
  vector<Trajectory> smoothed;
  if ( t.isValid() ) 
    { 
      vector<Trajectory> fitted = fitter()->fit(t);
      smoothingStep(fitted, smoothed);
    }
  return smoothed;
}
const TrajectoryFitter* KFFittingSmoother::fitter ( void  ) const [inline]

Definition at line 52 of file KFFittingSmoother.h.

References theFitter.

Referenced by KalmanAlignmentAlgorithm::initializeAlignmentSetups().

{return theFitter;}
const TrajectorySmoother* KFFittingSmoother::smoother ( ) const [inline]

Definition at line 53 of file KFFittingSmoother.h.

References theSmoother.

Referenced by KalmanAlignmentAlgorithm::initializeAlignmentSetups().

{return theSmoother;}
void KFFittingSmoother::smoothingStep ( std::vector< Trajectory > &  fitted,
std::vector< Trajectory > &  smoothed 
) const [private]

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

Definition at line 72 of file KFFittingSmoother.h.

Referenced by clone().

Definition at line 73 of file KFFittingSmoother.h.

Referenced by clone().

Definition at line 71 of file KFFittingSmoother.h.

Referenced by clone().

Definition at line 66 of file KFFittingSmoother.h.

Referenced by clone().

Definition at line 64 of file KFFittingSmoother.h.

Referenced by clone(), and fitter().

Definition at line 68 of file KFFittingSmoother.h.

Referenced by clone().

Definition at line 70 of file KFFittingSmoother.h.

Referenced by clone().

Definition at line 65 of file KFFittingSmoother.h.

Referenced by clone(), and smoother().

Definition at line 76 of file KFFittingSmoother.h.