CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/TrackingTools/GsfTracking/src/GsfTrajectorySmoother.cc

Go to the documentation of this file.
00001 #include "TrackingTools/GsfTracking/interface/GsfTrajectorySmoother.h"
00002 
00003 #include "TrackingTools/GsfTracking/interface/MultiTrajectoryStateMerger.h"
00004 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
00005 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
00006 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
00007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00008 
00009 GsfTrajectorySmoother::GsfTrajectorySmoother(const GsfPropagatorWithMaterial& aPropagator,
00010                                              const TrajectoryStateUpdator& aUpdator,
00011                                              const MeasurementEstimator& aEstimator,
00012                                              const MultiTrajectoryStateMerger& aMerger,
00013                                              float errorRescaling,
00014                                              const bool materialBeforeUpdate,
00015                                              const DetLayerGeometry* detLayerGeometry) :
00016   thePropagator(aPropagator.clone()),
00017   theGeomPropagator(0),
00018   theConvolutor(0),
00019   theUpdator(aUpdator.clone()),
00020   theEstimator(aEstimator.clone()),
00021   theMerger(aMerger.clone()),
00022   theMatBeforeUpdate(materialBeforeUpdate),
00023   theErrorRescaling(errorRescaling),
00024   theGeometry(detLayerGeometry)
00025 {
00026   if ( !theMatBeforeUpdate ) {
00027     theGeomPropagator = new GsfPropagatorAdapter(thePropagator->geometricalPropagator());
00028     theConvolutor = thePropagator->convolutionWithMaterial().clone();
00029   }
00030 
00031   if(!theGeometry) theGeometry = &dummyGeometry;
00032 
00033   //   static SimpleConfigurable<bool> timeConf(false,"GsfTrajectorySmoother:activateTiming");
00034   //   theTiming = timeConf.value();
00035 }
00036 
00037 GsfTrajectorySmoother::~GsfTrajectorySmoother() {
00038   delete thePropagator;
00039   delete theGeomPropagator;
00040   delete theConvolutor;
00041   delete theUpdator;
00042   delete theEstimator;
00043   delete theMerger;
00044 }
00045 
00046 std::vector<Trajectory> 
00047 GsfTrajectorySmoother::trajectories(const Trajectory& aTraj) const {
00048 
00049   //   static TimingReport::Item* propTimer =
00050   //     &(*TimingReport::current())[string("GsfTrajectorySmoother:propagation")];
00051   //   propTimer->switchCPU(false);
00052   //   if ( !theTiming )  propTimer->switchOn(false);
00053   //   static TimingReport::Item* updateTimer =
00054   //     &(*TimingReport::current())[string("GsfTrajectorySmoother:update")];
00055   //   updateTimer->switchCPU(false);
00056   //   if ( !theTiming )  updateTimer->switchOn(false);
00057   
00058   if(aTraj.empty()) return std::vector<Trajectory>();
00059   
00060   if (  aTraj.direction() == alongMomentum) {
00061     thePropagator->setPropagationDirection(oppositeToMomentum);
00062   }
00063   else {
00064     thePropagator->setPropagationDirection(alongMomentum);
00065   }
00066 
00067   Trajectory myTraj(aTraj.seed(), propagator()->propagationDirection());
00068   
00069   std::vector<TM> avtm = aTraj.measurements();
00070   
00071   TSOS predTsos = avtm.back().forwardPredictedState();
00072   predTsos.rescaleError(theErrorRescaling);
00073 
00074   if(!predTsos.isValid()) {
00075     edm::LogInfo("GsfTrajectorySmoother") 
00076       << "GsfTrajectorySmoother: predicted tsos of last measurement not valid!";
00077     return std::vector<Trajectory>();
00078   }
00079   TSOS currTsos;
00080 
00081   //first smoothed tm is last fitted
00082   if(avtm.back().recHit()->isValid()) {
00083     {
00084       //       TimeMe t(*updateTimer,false);
00085       currTsos = updator()->update(predTsos, *avtm.back().recHit());
00086     }
00087     if(!currTsos.isValid()) {
00088       edm::LogInfo("GsfTrajectorySmoother") << "GsfTrajectorySmoother: tsos not valid after update!";
00089       return std::vector<Trajectory>();
00090     }
00091 
00092     //check validity
00093     if (!avtm.back().forwardPredictedState().isValid() || !predTsos.isValid() || !avtm.back().updatedState().isValid()){
00094       edm::LogError("InvalidState")<<"first hit";
00095       return std::vector<Trajectory>();
00096     }
00097 
00098     myTraj.push(TM(avtm.back().forwardPredictedState(), 
00099                    predTsos,
00100                    avtm.back().updatedState(),
00101                    avtm.back().recHit(),
00102                    avtm.back().estimate(),
00103                    theGeometry->idToLayer(avtm.back().recHit()->geographicalId()) ), 
00104                 avtm.back().estimate());
00105   } else {
00106     currTsos = predTsos;
00107     //check validity
00108     if (!avtm.back().forwardPredictedState().isValid()){
00109       edm::LogError("InvalidState")<<"first hit on invalid hit";
00110       return std::vector<Trajectory>();
00111     }
00112 
00113     myTraj.push(TM(avtm.back().forwardPredictedState(),
00114                    avtm.back().recHit(),
00115                    0.,
00116                    theGeometry->idToLayer(avtm.back().recHit()->geographicalId() )  ));
00117   }
00118   
00119   TrajectoryStateCombiner combiner;
00120 
00121   for(std::vector<TM>::reverse_iterator itm = avtm.rbegin() + 1; 
00122       itm < avtm.rend() - 1; ++itm) {
00123     {
00124       //       TimeMe t(*propTimer,false);
00125       //       //
00126       //       // check type of surface in case of invalid hit
00127       //       // (in this version only propagations to planes are
00128       //       // supported for multi trajectory states)
00129       //       //
00130       //       if ( !(*itm).recHit().isValid() ) {
00131       //        const BoundPlane* plane = 
00132       //          dynamic_cast<const BoundPlane*>(&(*itm).recHit().det().surface());
00133       //        //
00134       //        // no plane: insert invalid 
00135       //        if ( plane==0 ) {
00136       //          myTraj.push(TM(TrajectoryStateOnSurface(),
00137       //                         (*itm).recHit());
00138       //          continue;
00139       //        }
00140       //     }
00141       predTsos = propagator()->propagate(currTsos,
00142                                          *(*itm).recHit()->surface());
00143     }
00144     if ( predTsos.isValid() && theConvolutor && theMatBeforeUpdate )
00145       predTsos = (*theConvolutor)(predTsos,
00146                                   propagator()->propagationDirection());
00147     if(!predTsos.isValid()) {
00148       edm::LogInfo("GsfTrajectorySmoother") << "GsfTrajectorySmoother: predicted tsos not valid!";
00149       return std::vector<Trajectory>();
00150     }
00151     if ( theMerger )  predTsos = theMerger->merge(predTsos);
00152 
00153     if((*itm).recHit()->isValid()) {
00154       //update
00155       {
00156         //      TimeMe t(*updateTimer,false);
00157         currTsos = updator()->update(predTsos, *(*itm).recHit());
00158       }
00159       if ( currTsos.isValid() && theConvolutor && !theMatBeforeUpdate )
00160         currTsos = (*theConvolutor)(currTsos,
00161                                     propagator()->propagationDirection());
00162       if(!currTsos.isValid()) {
00163         edm::LogInfo("GsfTrajectorySmoother") 
00164           << "GsfTrajectorySmoother: tsos not valid after update / material effects!";
00165         return std::vector<Trajectory>();
00166       }
00167       //3 different possibilities to calculate smoothed state:
00168       //1: update combined predictions with hit
00169       //2: combine fwd-prediction with bwd-filter
00170       //3: combine bwd-prediction with fwd-filter
00171       TSOS combTsos = combiner(predTsos, (*itm).forwardPredictedState());
00172       if(!combTsos.isValid()) {
00173         LogDebug("GsfTrajectorySmoother") << 
00174           "KFTrajectorySmoother: combined tsos not valid!\n"<<
00175           "pred Tsos pos: "<<predTsos.globalPosition()<< "\n" <<
00176           "pred Tsos mom: "<<predTsos.globalMomentum()<< "\n" <<
00177           "TrackingRecHit: "<<(*itm).recHit()->surface()->toGlobal((*itm).recHit()->localPosition())<< "\n" ;
00178         return std::vector<Trajectory>();
00179       }
00180 
00181       TSOS smooTsos = combiner((*itm).updatedState(), predTsos);
00182 
00183       if(!smooTsos.isValid()) {
00184         LogDebug("GsfTrajectorySmoother") <<
00185           "KFTrajectorySmoother: smoothed tsos not valid!";
00186         return std::vector<Trajectory>();
00187       }
00188 
00189       if (!(*itm).forwardPredictedState().isValid() || !predTsos.isValid() || !smooTsos.isValid() ){
00190         edm::LogError("InvalidState")<<"inside hits with combination.";
00191         return std::vector<Trajectory>();
00192       }
00193 
00194 
00195       myTraj.push(TM((*itm).forwardPredictedState(),
00196                      predTsos,
00197                      smooTsos,
00198                      (*itm).recHit(),
00199                      estimator()->estimate(combTsos, *(*itm).recHit()).second,
00200                      theGeometry->idToLayer((*itm).recHit()->geographicalId() ) ),
00201                   (*itm).estimate());
00202     } 
00203     else {
00204       currTsos = predTsos;
00205       TSOS combTsos = combiner(predTsos, (*itm).forwardPredictedState());
00206       
00207       if(!combTsos.isValid()) {
00208         LogDebug("GsfTrajectorySmoother") << 
00209           "KFTrajectorySmoother: combined tsos not valid!";
00210         return std::vector<Trajectory>();
00211       }
00212 
00213       if (!(*itm).forwardPredictedState().isValid() || !predTsos.isValid() || !combTsos.isValid() ){
00214         edm::LogError("InvalidState")<<"inside hits with invalid rechit.";
00215         return std::vector<Trajectory>();
00216       }
00217 
00218       myTraj.push(TM((*itm).forwardPredictedState(),
00219                      predTsos,
00220                      combTsos,
00221                      (*itm).recHit(),
00222                      0.,
00223                      theGeometry->idToLayer((*itm).recHit()->geographicalId()) ));
00224     }
00225     if ( theMerger )  currTsos = theMerger->merge(currTsos);
00226   }
00227 
00228   //last smoothed tm is last filtered
00229   {
00230     //     TimeMe t(*propTimer,false);
00231     predTsos = propagator()->propagate(currTsos,
00232                                        *avtm.front().recHit()->surface());
00233   }
00234   if ( predTsos.isValid() && theConvolutor && theMatBeforeUpdate )
00235     predTsos = (*theConvolutor)(predTsos,
00236                                 propagator()->propagationDirection());
00237   if(!predTsos.isValid()) {
00238     edm::LogInfo("GsfTrajectorySmoother") << "GsfTrajectorySmoother: predicted tsos not valid!";
00239     return std::vector<Trajectory>();
00240   }
00241   if ( theMerger )  predTsos = theMerger->merge(predTsos);
00242 
00243   if(avtm.front().recHit()->isValid()) {
00244     //update
00245     {
00246       //       TimeMe t(*updateTimer,false);
00247       currTsos = updator()->update(predTsos, *avtm.front().recHit());
00248     }
00249     if ( currTsos.isValid() && theConvolutor && !theMatBeforeUpdate )
00250       currTsos = (*theConvolutor)(currTsos,
00251                                   propagator()->propagationDirection());
00252     if(!currTsos.isValid()) {
00253       edm::LogInfo("GsfTrajectorySmoother") 
00254         << "GsfTrajectorySmoother: tsos not valid after update / material effects!";
00255       return std::vector<Trajectory>();
00256     }
00257   
00258     if (!avtm.front().forwardPredictedState().isValid() || !predTsos.isValid() || !currTsos.isValid() ){
00259       edm::LogError("InvalidState")<<"last hit";
00260       return std::vector<Trajectory>();
00261     }
00262 
00263     myTraj.push(TM(avtm.front().forwardPredictedState(),
00264                    predTsos,
00265                    currTsos,
00266                    avtm.front().recHit(),
00267                    estimator()->estimate(predTsos, *avtm.front().recHit()).second,
00268                    theGeometry->idToLayer(avtm.front().recHit()->geographicalId() )),
00269                 avtm.front().estimate());
00270     //estimator()->estimate(predTsos, avtm.front().recHit()));
00271   }
00272   else {
00273     if (!avtm.front().forwardPredictedState().isValid()){
00274       edm::LogError("InvalidState")<<"last invalid hit";
00275       return std::vector<Trajectory>();
00276     }
00277     myTraj.push(TM(avtm.front().forwardPredictedState(),
00278                    avtm.front().recHit(),
00279                    0.,
00280                    theGeometry->idToLayer(avtm.front().recHit()->geographicalId()) ));
00281   }
00282 
00283   return std::vector<Trajectory>(1, myTraj); 
00284 }