CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoMuon/TrackingTools/src/MuonTrajectoryUpdator.cc

Go to the documentation of this file.
00001 
00017 #include "RecoMuon/TrackingTools/interface/MuonTrajectoryUpdator.h"
00018 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
00019 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00020 
00021 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
00022 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
00023 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00025 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
00026 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
00027 #include "TrackingTools/DetLayers/interface/DetLayer.h"
00028 
00029 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
00030 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBreaker.h"
00031 
00032 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00033 #include "FWCore/Utilities/interface/Exception.h"
00034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00035 
00036 #include <algorithm>
00037 
00038 using namespace edm;
00039 using namespace std;
00040 
00042 MuonTrajectoryUpdator::MuonTrajectoryUpdator(const edm::ParameterSet& par,
00043                                              NavigationDirection fitDirection): theFitDirection(fitDirection){
00044   
00045   // The max allowed chi2 to accept a rechit in the fit
00046   theMaxChi2 = par.getParameter<double>("MaxChi2");
00047   theEstimator = new Chi2MeasurementEstimator(theMaxChi2);
00048   
00049   // The KF updator
00050   theUpdator= new KFUpdator();
00051 
00052   // The granularity
00053   theGranularity = par.getParameter<int>("Granularity");
00054 
00055   // Rescale the error of the first state?
00056   theRescaleErrorFlag =  par.getParameter<bool>("RescaleError");
00057 
00058   if(theRescaleErrorFlag)
00059     // The rescale factor
00060     theRescaleFactor =  par.getParameter<double>("RescaleErrorFactor");
00061   
00062   // Use the invalid hits?
00063   useInvalidHits =  par.getParameter<bool>("UseInvalidHits");
00064 
00065   // Flag needed for the rescaling
00066   theFirstTSOSFlag = true;
00067 
00068   // Exlude RPC from the fit?
00069    theRPCExFlag = par.getParameter<bool>("ExcludeRPCFromFit");
00070 }
00071 
00072 MuonTrajectoryUpdator::MuonTrajectoryUpdator( NavigationDirection fitDirection,
00073                                               double chi2, int granularity): theMaxChi2(chi2),
00074                                                                              theGranularity(granularity),
00075                                                                              theFitDirection(fitDirection){
00076   theEstimator = new Chi2MeasurementEstimator(theMaxChi2);
00077   
00078   // The KF updator
00079   theUpdator= new KFUpdator();
00080 }
00081 
00083 MuonTrajectoryUpdator::~MuonTrajectoryUpdator(){
00084   delete theEstimator;
00085   delete theUpdator;
00086 }
00087 
00088 void MuonTrajectoryUpdator::makeFirstTime(){
00089   theFirstTSOSFlag = true;
00090 }
00091 
00092 
00093 pair<bool,TrajectoryStateOnSurface> 
00094 MuonTrajectoryUpdator::update(const TrajectoryMeasurement* measurement,
00095                               Trajectory& trajectory,
00096                               const Propagator *propagator){
00097   
00098   const std::string metname = "Muon|RecoMuon|MuonTrajectoryUpdator";
00099 
00100   MuonPatternRecoDumper muonDumper;
00101 
00102   // Status of the updating
00103   bool updated=false;
00104   
00105   if(!measurement) return pair<bool,TrajectoryStateOnSurface>(updated,TrajectoryStateOnSurface() );
00106 
00107   // measurement layer
00108   const DetLayer* detLayer=measurement->layer();
00109 
00110   // these are the 4D segment for the CSC/DT and a point for the RPC 
00111   TransientTrackingRecHit::ConstRecHitPointer muonRecHit = measurement->recHit();
00112  
00113   // The KFUpdator takes TransientTrackingRecHits as arg.
00114   TransientTrackingRecHit::ConstRecHitContainer recHitsForFit =
00115   MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit,theGranularity);
00116 
00117   // sort the container in agreement with the porpagation direction
00118   sort(recHitsForFit,detLayer);
00119   
00120   TrajectoryStateOnSurface lastUpdatedTSOS = measurement->predictedState();
00121   
00122   LogTrace(metname)<<"Number of rechits for the fit: "<<recHitsForFit.size()<<endl;
00123  
00124   TransientTrackingRecHit::ConstRecHitContainer::iterator recHit;
00125   for(recHit = recHitsForFit.begin(); recHit != recHitsForFit.end(); ++recHit ) {
00126     if ((*recHit)->isValid() ) {
00127 
00128       // propagate the TSOS onto the rechit plane
00129       TrajectoryStateOnSurface propagatedTSOS  = propagateState(lastUpdatedTSOS, measurement, 
00130                                                                 *recHit, propagator);
00131       
00132       if ( propagatedTSOS.isValid() ) {
00133         pair<bool,double> thisChi2 = estimator()->estimate(propagatedTSOS, *((*recHit).get()));
00134 
00135         LogTrace(metname) << "Estimation for Kalman Fit. Chi2: " << thisChi2.second;
00136 
00137         // Is an RPC hit? Prepare the variable to possibly exluding it from the fit
00138         bool wantIncludeThisHit = true;
00139         if (theRPCExFlag && 
00140             (*recHit)->geographicalId().det() == DetId::Muon &&
00141             (*recHit)->geographicalId().subdetId() == MuonSubdetId::RPC){
00142           wantIncludeThisHit = false;   
00143           LogTrace(metname) << "This is an RPC hit and the present configuration is such that it will be excluded from the fit";
00144         }
00145 
00146         
00147         // The Chi2 cut was already applied in the estimator, which
00148         // returns 0 if the chi2 is bigger than the cut defined in its
00149         // constructor
00150         if (thisChi2.first) {
00151           updated=true;
00152           if (wantIncludeThisHit) { // This split is a trick to have the RPC hits counted as updatable (in used chamber counting), while are not actually included in the fit when the proper obtion is activated.
00153 
00154           LogTrace(metname) << endl 
00155                             << "     Kalman Start" << "\n" << "\n";
00156           LogTrace(metname) << "  Meas. Position : " << (**recHit).globalPosition() << "\n"
00157                             << "  Pred. Position : " << propagatedTSOS.globalPosition()
00158                             << "  Pred Direction : " << propagatedTSOS.globalDirection()<< endl;
00159 
00160           if(theRescaleErrorFlag && theFirstTSOSFlag){
00161             propagatedTSOS.rescaleError(theRescaleFactor);
00162             theFirstTSOSFlag = false;
00163           }
00164 
00165           lastUpdatedTSOS = measurementUpdator()->update(propagatedTSOS,*((*recHit).get()));
00166 
00167           LogTrace(metname) << "  Fit   Position : " << lastUpdatedTSOS.globalPosition()
00168                             << "  Fit  Direction : " << lastUpdatedTSOS.globalDirection()
00169                             << "\n"
00170                             << "  Fit position radius : " 
00171                             << lastUpdatedTSOS.globalPosition().perp()
00172                             << "filter updated" << endl;
00173           
00174           LogTrace(metname) << muonDumper.dumpTSOS(lastUpdatedTSOS);
00175           
00176           LogTrace(metname) << "\n\n     Kalman End" << "\n" << "\n";         
00177           
00178           TrajectoryMeasurement updatedMeasurement = updateMeasurement( propagatedTSOS, lastUpdatedTSOS, 
00179                                                                         *recHit, thisChi2.second, detLayer, 
00180                                                                         measurement);
00181           // FIXME: check!
00182           trajectory.push(updatedMeasurement, thisChi2.second); 
00183           }
00184           else {
00185             LogTrace(metname) << "  Compatible RecHit with good chi2 but made with RPC when it was decided to not include it in the fit"
00186                               << "  --> trajectory NOT updated, invalid RecHit added." << endl;
00187               
00188             MuonTransientTrackingRecHit::MuonRecHitPointer invalidRhPtr = MuonTransientTrackingRecHit::specificBuild( (*recHit)->det(), (*recHit)->hit() );
00189             invalidRhPtr->invalidateHit();
00190             TrajectoryMeasurement invalidRhMeasurement(propagatedTSOS, propagatedTSOS, invalidRhPtr.get(), thisChi2.second, detLayer);
00191             trajectory.push(invalidRhMeasurement, thisChi2.second);                 
00192           }
00193         }
00194         else {
00195           if(useInvalidHits) {
00196             LogTrace(metname) << "  Compatible RecHit with too large chi2"
00197                             << "  --> trajectory NOT updated, invalid RecHit added." << endl;
00198 
00199             MuonTransientTrackingRecHit::MuonRecHitPointer invalidRhPtr = MuonTransientTrackingRecHit::specificBuild( (*recHit)->det(), (*recHit)->hit() );
00200             invalidRhPtr->invalidateHit();
00201             TrajectoryMeasurement invalidRhMeasurement(propagatedTSOS, propagatedTSOS, invalidRhPtr.get(), thisChi2.second, detLayer);
00202             trajectory.push(invalidRhMeasurement, thisChi2.second);       
00203           }
00204         }
00205       }
00206     }
00207   }
00208   recHitsForFit.clear();
00209   return pair<bool,TrajectoryStateOnSurface>(updated,lastUpdatedTSOS);
00210 }
00211 
00212 TrajectoryStateOnSurface 
00213 MuonTrajectoryUpdator::propagateState(const TrajectoryStateOnSurface& state,
00214                                       const TrajectoryMeasurement* measurement, 
00215                                       const TransientTrackingRecHit::ConstRecHitPointer  & current,
00216                                       const Propagator *propagator) const{
00217 
00218   const TransientTrackingRecHit::ConstRecHitPointer mother = measurement->recHit();
00219 
00220   if( current->geographicalId() == mother->geographicalId() )
00221     return measurement->predictedState();
00222   
00223   const TrajectoryStateOnSurface  tsos =
00224     propagator->propagate(state, current->det()->surface());
00225   return tsos;
00226 
00227 }
00228 
00229 // FIXME: would I a different threatment for the two prop dirrections??
00230 TrajectoryMeasurement MuonTrajectoryUpdator::updateMeasurement(  const TrajectoryStateOnSurface &propagatedTSOS, 
00231                                                                  const TrajectoryStateOnSurface &lastUpdatedTSOS, 
00232                                                                  const TransientTrackingRecHit::ConstRecHitPointer &recHit,
00233                                                                  const double &chi2, const DetLayer *detLayer, 
00234                                                                  const TrajectoryMeasurement *initialMeasurement){
00235   return TrajectoryMeasurement(propagatedTSOS, lastUpdatedTSOS, 
00236                                recHit,chi2,detLayer);
00237 
00238   //   // FIXME: put a better check! One could fit in first out-in and then in - out 
00239   //   if(propagator()->propagationDirection() == alongMomentum) 
00240   //     return TrajectoryMeasurement(propagatedTSOS, lastUpdatedTSOS, 
00241   //                             recHit,thisChi2.second,detLayer);
00242   
00243   //   // FIXME: Check this carefully!!
00244   //   else if(propagator()->propagationDirection() == oppositeToMomentum)
00245   //     return TrajectoryMeasurement(initialMeasurement->forwardPredictedState(),
00246   //                             propagatedTSOS, lastUpdatedTSOS, 
00247   //                             recHit,thisChi2.second,detLayer);
00248   //   else{
00249   //     LogError("MuonTrajectoryUpdator::updateMeasurement") <<"Wrong propagation direction!!";
00250   //   }
00251 }
00252 
00253 
00254 void MuonTrajectoryUpdator::sort(TransientTrackingRecHit::ConstRecHitContainer& recHitsForFit, 
00255                                  const DetLayer* detLayer){
00256   
00257   if(detLayer->subDetector()==GeomDetEnumerators::DT){
00258     if(fitDirection() == insideOut)
00259       stable_sort(recHitsForFit.begin(),recHitsForFit.end(), RadiusComparatorInOut() );
00260     else if(fitDirection() == outsideIn)
00261       stable_sort(recHitsForFit.begin(),recHitsForFit.end(),RadiusComparatorOutIn() ); 
00262     else
00263       LogError("Muon|RecoMuon|MuonTrajectoryUpdator") <<"MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
00264   }
00265 
00266   else if(detLayer->subDetector()==GeomDetEnumerators::CSC){
00267     if(fitDirection() == insideOut)
00268       stable_sort(recHitsForFit.begin(),recHitsForFit.end(), ZedComparatorInOut() );
00269     else if(fitDirection() == outsideIn)
00270       stable_sort(recHitsForFit.begin(),recHitsForFit.end(), ZedComparatorOutIn() );  
00271     else
00272       LogError("Muon|RecoMuon|MuonTrajectoryUpdator") <<"MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
00273   }
00274 }