CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc

Go to the documentation of this file.
00001 #include "TrackingTools/GsfTracking/interface/GsfTrajectoryFitter.h"
00002 
00003 #include "TrackingTools/TrackFitters/interface/TrajectoryStateWithArbitraryError.h"
00004 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
00005 // #include "CommonDet/BasicDet/interface/Det.h"
00006 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
00007 #include "TrackingTools/GsfTracking/interface/MultiTrajectoryStateMerger.h"
00008 // #include "Utilities/Notification/interface/Verbose.h"
00009 // #include "Utilities/Notification/interface/TimingReport.h"
00010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00011 
00012 GsfTrajectoryFitter::GsfTrajectoryFitter(const Propagator& aPropagator,
00013                                          const TrajectoryStateUpdator& aUpdator,
00014                                          const MeasurementEstimator& aEstimator,
00015                                          const MultiTrajectoryStateMerger& aMerger,
00016                                          const DetLayerGeometry* detLayerGeometry) :
00017   thePropagator(aPropagator.clone()),
00018   theUpdator(aUpdator.clone()),
00019   theEstimator(aEstimator.clone()),
00020   theMerger(aMerger.clone()),
00021   theGeometry(detLayerGeometry)
00022 {
00023   if(!theGeometry) theGeometry = &dummyGeometry;
00024   //   static SimpleConfigurable<bool> timeConf(false,"GsfTrajectoryFitter:activateTiming");
00025   //   theTiming = timeConf.value();
00026 }
00027 
00028 GsfTrajectoryFitter::~GsfTrajectoryFitter() {
00029   delete thePropagator;
00030   delete theUpdator;
00031   delete theEstimator;
00032   delete theMerger;
00033 }
00034 
00035 std::vector<Trajectory> GsfTrajectoryFitter::fit(const Trajectory& aTraj) const 
00036 {  
00037   if(aTraj.empty()) return std::vector<Trajectory>();
00038  
00039   TM firstTM = aTraj.firstMeasurement();
00040   TSOS firstTsos = TrajectoryStateWithArbitraryError()(firstTM.updatedState());
00041   
00042   return fit(aTraj.seed(), aTraj.recHits(), firstTsos);
00043 }
00044 
00045 std::vector<Trajectory> GsfTrajectoryFitter::fit(const TrajectorySeed& aSeed,
00046                                                  const RecHitContainer& hits) const {
00047 
00048   edm::LogError("GsfTrajectoryFitter") 
00049     << "GsfTrajectoryFitter::fit(TrajectorySeed, vector<RecHit>) not implemented";
00050 
00051   return std::vector<Trajectory>();
00052 }
00053 
00054 std::vector<Trajectory> GsfTrajectoryFitter::fit(const TrajectorySeed& aSeed,
00055                                                  const RecHitContainer& hits, 
00056                                                  const TSOS& firstPredTsos) const {
00057 
00058   //   static TimingReport::Item* propTimer =
00059   //     &(*TimingReport::current())[string("GsfTrajectoryFitter:propagation")];
00060   //   propTimer->switchCPU(false);
00061   //   if ( !theTiming )  propTimer->switchOn(false);
00062   //   static TimingReport::Item* updateTimer =
00063   //     &(*TimingReport::current())[string("GsfTrajectoryFitter:update")];
00064   //   updateTimer->switchCPU(false);
00065   //   if ( !theTiming )  updateTimer->switchOn(false);
00066 
00067   if(hits.empty()) return std::vector<Trajectory>();
00068 
00069   Trajectory myTraj(aSeed, propagator()->propagationDirection());
00070 
00071   TSOS predTsos(firstPredTsos);
00072   if(!predTsos.isValid()) {
00073     edm::LogInfo("GsfTrajectoryFitter") 
00074       << "GsfTrajectoryFitter: predicted tsos of first measurement not valid!";
00075     return std::vector<Trajectory>();
00076   } 
00077 
00078   TSOS currTsos;
00079   if(hits.front()->isValid()) {
00080     //update
00081     TransientTrackingRecHit::RecHitPointer preciseHit = hits.front()->clone(predTsos);
00082     {
00083       //       TimeMe t(*updateTimer,false);
00084       currTsos = updator()->update(predTsos, *preciseHit);
00085     }
00086     if (!predTsos.isValid() || !currTsos.isValid()){
00087       edm::LogError("InvalidState")<<"first hit";
00088       return std::vector<Trajectory>();
00089     }
00090     myTraj.push(TM(predTsos, currTsos, preciseHit, 0., theGeometry->idToLayer(preciseHit->geographicalId() )),
00091                 estimator()->estimate(predTsos, *preciseHit).second);
00092   } else {
00093     currTsos = predTsos;
00094     if (!predTsos.isValid()){
00095       edm::LogError("InvalidState")<<"first invalid hit";
00096       return std::vector<Trajectory>();
00097     }
00098     myTraj.push(TM(predTsos, *hits.begin(),0., theGeometry->idToLayer((*hits.begin())->geographicalId()) ));
00099   }
00100   
00101   for(RecHitContainer::const_iterator ihit = hits.begin() + 1; 
00102       ihit != hits.end(); ihit++) {    
00103     //
00104     // temporary protection copied from KFTrajectoryFitter.
00105     //
00106     if ((**ihit).isValid() == false && (**ihit).det() == 0) {
00107       LogDebug("GsfTrajectoryFitter") << " Error: invalid hit with no GeomDet attached .... skipping";
00108       continue;
00109     }
00110 
00112     //     //
00113     //     // check type of surface in case of invalid hit
00114     //     // (in this version only propagations to planes are
00115     //     // supported for multi trajectory states)
00116     //     //
00117     //     if ( !(**ihit).isValid() ) {
00118     //       const BoundPlane* plane = 
00119     //  dynamic_cast<const BoundPlane*>(&(**ihit).det().surface());
00120     //       //
00121     //       // no plane: insert invalid measurement
00122     //       //
00123     //       if ( plane==0 ) {
00124     //  myTraj.push(TM(TrajectoryStateOnSurface(),&(**ihit)));
00125     //  continue;
00126     //       }
00127     //     }
00128     {
00129       //       TimeMe t(*propTimer,false);
00130       predTsos = propagator()->propagate(currTsos,
00131                                          (**ihit).det()->surface());
00132     }
00133     if(!predTsos.isValid()) {
00134       if ( myTraj.foundHits()>=3 ) {
00135         edm::LogInfo("GsfTrajectoryFitter") 
00136           << "GsfTrajectoryFitter: predicted tsos not valid! \n"
00137           << "Returning trajectory with " << myTraj.foundHits() << " found hits.";
00138         return std::vector<Trajectory>(1,myTraj);
00139       }
00140       else {
00141       edm::LogInfo("GsfTrajectoryFitter") 
00142         << "GsfTrajectoryFitter: predicted tsos not valid after " << myTraj.foundHits()
00143         << " hits, discarding candidate!";
00144         return std::vector<Trajectory>();
00145       }
00146     }
00147     if ( merger() ) predTsos = merger()->merge(predTsos);
00148     
00149     if((**ihit).isValid()) {
00150       //update
00151       TransientTrackingRecHit::RecHitPointer preciseHit = (**ihit).clone(predTsos);
00152       {
00153         //      TimeMe t(*updateTimer,false);
00154         currTsos = updator()->update(predTsos, *preciseHit);
00155       }
00156       if (!predTsos.isValid() || !currTsos.isValid()){
00157         edm::LogError("InvalidState")<<"inside hit";
00158         return std::vector<Trajectory>();
00159       }
00160       myTraj.push(TM(predTsos, currTsos, preciseHit,
00161                      estimator()->estimate(predTsos, *preciseHit).second,
00162                      theGeometry->idToLayer(preciseHit->geographicalId() )));
00163     } else {
00164       currTsos = predTsos;
00165       if (!predTsos.isValid()){
00166       edm::LogError("InvalidState")<<"inside invalid hit";
00167       return std::vector<Trajectory>();
00168       }
00169       myTraj.push(TM(predTsos, *ihit,0., theGeometry->idToLayer( (*ihit)->geographicalId()) ));
00170     }
00171   }
00172   return std::vector<Trajectory>(1, myTraj);
00173 }