CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/TrackingTools/TrackRefitter/src/TrackTransformerForCosmicMuons.cc

Go to the documentation of this file.
00001 #include "TrackingTools/TrackRefitter/interface/TrackTransformerForCosmicMuons.h"
00002 
00003 #include "FWCore/Framework/interface/EventSetup.h"
00004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 
00007 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00008 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00009 
00010 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00011 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
00012 #include "TrackingTools/PatternTools/interface/TrajectorySmoother.h"
00013 
00014 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00015 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00016 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
00017 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
00018 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00019 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
00020 
00021 #include "DataFormats/TrackReco/interface/Track.h"
00022 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00023 #include "DataFormats/DetId/interface/DetId.h"
00024 
00025 
00026 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
00027 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
00028 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
00029 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
00030 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00031 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00032 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
00033 #include "DataFormats/MuonDetId/interface/DTWireId.h"
00034 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
00035 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00036 
00037 
00038 using namespace std;
00039 using namespace edm;
00040 
00042 TrackTransformerForCosmicMuons::TrackTransformerForCosmicMuons(const ParameterSet& parameterSet){
00043   
00044   theTrackerRecHitBuilderName = parameterSet.getParameter<string>("TrackerRecHitBuilder");
00045   theMuonRecHitBuilderName = parameterSet.getParameter<string>("MuonRecHitBuilder");
00046 
00047   theRPCInTheFit = parameterSet.getParameter<bool>("RefitRPCHits");
00048 
00049   theCacheId_TC = theCacheId_GTG = theCacheId_MG = theCacheId_TRH = 0;
00050 }
00051 
00053 TrackTransformerForCosmicMuons::~TrackTransformerForCosmicMuons(){}
00054 
00055 
00056 void TrackTransformerForCosmicMuons::setServices(const EventSetup& setup){
00057   
00058   const std::string metname = "Reco|TrackingTools|TrackTransformer";
00059 
00060   setup.get<TrajectoryFitter::Record>().get("KFFitterForRefitInsideOut",theFitterIO);
00061   setup.get<TrajectoryFitter::Record>().get("KFSmootherForRefitInsideOut",theSmootherIO);  
00062   setup.get<TrajectoryFitter::Record>().get("KFFitterForRefitOutsideIn",theFitterOI);
00063   setup.get<TrajectoryFitter::Record>().get("KFSmootherForRefitOutsideIn",theSmootherOI);
00064 
00065   unsigned long long newCacheId_TC = setup.get<TrackingComponentsRecord>().cacheIdentifier();
00066 
00067   if ( newCacheId_TC != theCacheId_TC ){
00068     LogTrace(metname) << "Tracking Component changed!";
00069     theCacheId_TC = newCacheId_TC;
00070     setup.get<TrackingComponentsRecord>().get("SmartPropagatorRK",thePropagatorIO);
00071     setup.get<TrackingComponentsRecord>().get("SmartPropagatorRKOpposite",thePropagatorOI);
00072   }
00073 
00074   // Global Tracking Geometry
00075   unsigned long long newCacheId_GTG = setup.get<GlobalTrackingGeometryRecord>().cacheIdentifier();
00076   if ( newCacheId_GTG != theCacheId_GTG ) {
00077     LogTrace(metname) << "GlobalTrackingGeometry changed!";
00078     theCacheId_GTG = newCacheId_GTG;
00079     setup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry); 
00080   }
00081   
00082   // Magfield Field
00083   unsigned long long newCacheId_MG = setup.get<IdealMagneticFieldRecord>().cacheIdentifier();
00084   if ( newCacheId_MG != theCacheId_MG ) {
00085     LogTrace(metname) << "Magnetic Field changed!";
00086     theCacheId_MG = newCacheId_MG;
00087     setup.get<IdealMagneticFieldRecord>().get(theMGField);
00088   }
00089   
00090   // Transient Rechit Builders
00091   unsigned long long newCacheId_TRH = setup.get<TransientRecHitRecord>().cacheIdentifier();
00092   if ( newCacheId_TRH != theCacheId_TRH ) {
00093     theCacheId_TRH = newCacheId_TRH;
00094     LogTrace(metname) << "TransientRecHitRecord changed!";
00095     setup.get<TransientRecHitRecord>().get(theTrackerRecHitBuilderName,theTrackerRecHitBuilder);
00096     setup.get<TransientRecHitRecord>().get(theMuonRecHitBuilderName,theMuonRecHitBuilder);
00097   }
00098 }
00099 
00100 
00101 TransientTrackingRecHit::ConstRecHitContainer
00102 TrackTransformerForCosmicMuons::getTransientRecHits(const reco::TransientTrack& track) const {
00103 
00104   TransientTrackingRecHit::ConstRecHitContainer tkHits;
00105   TransientTrackingRecHit::ConstRecHitContainer staHits;
00106 
00107   bool printout = false;
00108   
00109   bool quad1 = false;
00110   bool quad2 = false;
00111   bool quad3 = false;
00112   bool quad4 = false;
00113 
00114   for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit)
00115     if((*hit)->isValid())
00116       if ( (*hit)->geographicalId().det() == DetId::Muon ){
00117                 if( (*hit)->geographicalId().subdetId() == 3 && !theRPCInTheFit){
00118                         LogTrace("Reco|TrackingTools|TrackTransformer") << "RPC Rec Hit discarged"; 
00119                         continue;
00120                 }
00121                 staHits.push_back(theMuonRecHitBuilder->build(&**hit));
00122         DetId hitId = staHits.back()->geographicalId();
00123         GlobalPoint glbpoint = trackingGeometry()->idToDet(hitId)->position();
00124                 float gpy=glbpoint.y(); 
00125                 float gpz=glbpoint.z(); 
00126 //              if (gpy != 0 && gpz !=0) slopeSum += gpy / gpz;
00127 
00128                 if (gpy > 0 && gpz > 0) quad1 = true; 
00129                 else if (gpy > 0 && gpz < 0) quad2 = true; 
00130                 else if (gpy < 0 && gpz < 0) quad3 = true; 
00131                 else if (gpy < 0 && gpz > 0) quad4 = true; 
00132                 else return tkHits;
00133       }
00134 
00135   if(staHits.empty()) return staHits;
00136 
00137   if (quad1 && quad2) return tkHits;
00138   if (quad1 && quad3) return tkHits;
00139   if (quad1 && quad4) return tkHits;
00140   if (quad2 && quad3) return tkHits;
00141   if (quad2 && quad4) return tkHits;
00142   if (quad3 && quad4) return tkHits;
00143 
00144 
00145   bool doReverse = staHits.front()->globalPosition().y()>0 ? true : false;
00146 
00147   if (quad1) doReverse = SlopeSum(staHits); 
00148   if (quad2) doReverse = SlopeSum(staHits); 
00149   if (quad3) doReverse = SlopeSum(staHits); 
00150   if (quad4) doReverse = SlopeSum(staHits); 
00151   if(doReverse){
00152     reverse(staHits.begin(),staHits.end());
00153   }
00154 
00155   copy(staHits.begin(),staHits.end(),back_inserter(tkHits));
00156 
00157 
00159 //  if ( quad2 && slopeSum > 0) printout = true;
00162 //  swap = printout;
00163 
00164   printout = quad1;  
00165 
00166   if (printout) for(TransientTrackingRecHit::ConstRecHitContainer::const_iterator hit = tkHits.begin();
00167       hit !=tkHits.end(); ++hit){
00168 
00169     DetId hitId = (*hit)->geographicalId();
00170     GlobalPoint glbpoint = trackingGeometry()->idToDet(hitId)->position();
00171 
00172 
00173     if(hitId.det() == DetId::Muon) {
00174       if(hitId.subdetId() == MuonSubdetId::DT)
00175                 LogTrace("TrackFitters") << glbpoint << " I am DT " << DTWireId(hitId);
00176 //              std::cout<< glbpoint << " I am DT " << DTWireId(hitId)<<std::endl;
00177       else if (hitId.subdetId() == MuonSubdetId::CSC )
00178                 LogTrace("TrackFitters") << glbpoint << " I am CSC " << CSCDetId(hitId);
00179 //              std::cout<< glbpoint << " I am CSC " << CSCDetId(hitId)<<std::endl;
00180       else if (hitId.subdetId() == MuonSubdetId::RPC )
00181                 LogTrace("TrackFitters") << glbpoint << " I am RPC " << RPCDetId(hitId);
00182       else 
00183                 LogTrace("TrackFitters") << " UNKNOWN MUON HIT TYPE ";
00184     } 
00185   } 
00186   return tkHits;
00187 }
00188 
00189 
00191 ESHandle<TrajectoryFitter> TrackTransformerForCosmicMuons::fitter(bool up, int quad, float sumy) const{
00192   if(quad ==1) {if (sumy < 0) return theFitterOI; else return theFitterIO;}
00193   if(quad ==2) {if (sumy < 0) return theFitterOI; else return theFitterIO;}
00194   if(quad ==3) {if (sumy > 0) return theFitterOI; else return theFitterIO;}
00195   if(quad ==4) {if (sumy > 0) return theFitterOI; else return theFitterIO;}
00196 
00197   if(up) return theFitterOI;
00198   else return theFitterIO;
00199 }
00200   
00202 ESHandle<TrajectorySmoother> TrackTransformerForCosmicMuons::smoother(bool up, int quad, float sumy) const{
00203   if(quad ==1){ if (sumy < 0) return theSmootherOI; else return theSmootherIO;}
00204   if(quad ==2){ if (sumy < 0) return theSmootherOI; else return theSmootherIO;}
00205   if(quad ==3){ if (sumy > 0) return theSmootherOI; else return theSmootherIO;}
00206   if(quad ==4){ if (sumy > 0) return theSmootherOI; else return theSmootherIO;}
00207   if(up) return theSmootherOI;
00208   else return theSmootherIO;
00209 }
00210 
00211 ESHandle<Propagator> TrackTransformerForCosmicMuons::propagator(bool up, int quad, float sumy) const{
00212   if(quad ==1) {if (sumy > 0) return thePropagatorOI; else return thePropagatorIO;}
00213   if(quad ==2) {if (sumy > 0) return thePropagatorOI; else return thePropagatorIO;}
00214   if(quad ==3) {if (sumy < 0) return thePropagatorOI; else return thePropagatorIO;}
00215   if(quad ==4) {if (sumy < 0) return thePropagatorOI; else return thePropagatorIO;}
00216   if(up) return thePropagatorIO;
00217   else return thePropagatorOI;
00218 }
00219 
00220 
00221 
00223 vector<Trajectory> TrackTransformerForCosmicMuons::transform(const reco::Track& tr) const {
00224 
00225   const std::string metname = "Reco|TrackingTools|TrackTransformer";
00226   
00227   reco::TransientTrack track(tr,magneticField(),trackingGeometry());   
00228 
00229   // Build the transient Rechits
00230   TransientTrackingRecHit::ConstRecHitContainer recHitsForReFit;// = getTransientRecHits(track);
00231   TransientTrackingRecHit::ConstRecHitContainer staHits;
00232 
00233 
00234   bool quad1 = false;
00235   bool quad2 = false;
00236   bool quad3 = false;
00237   bool quad4 = false;
00238   int quadrant = 0;
00239 
00240   for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit)
00241     if((*hit)->isValid())
00242       if ( (*hit)->geographicalId().det() == DetId::Muon ){
00243                 if( (*hit)->geographicalId().subdetId() == 3 && !theRPCInTheFit){
00244                         LogTrace("Reco|TrackingTools|TrackTransformer") << "RPC Rec Hit discarged"; 
00245                         continue;
00246                 }
00247                 staHits.push_back(theMuonRecHitBuilder->build(&**hit));
00248         DetId hitId = staHits.back()->geographicalId();
00249         GlobalPoint glbpoint = trackingGeometry()->idToDet(hitId)->position();
00250                 float gpy=glbpoint.y(); 
00251                 float gpz=glbpoint.z(); 
00252 //              if (gpy != 0 && gpz !=0) slopeSum += gpy / gpz;
00253 
00254                 if (gpy > 0 && gpz > 0)         {quad1 = true; quadrant = 1;}
00255                 else if (gpy > 0 && gpz < 0){quad2 = true; quadrant = 2;}
00256                 else if (gpy < 0 && gpz < 0){quad3 = true; quadrant = 3;}
00257                 else if (gpy < 0 && gpz > 0){quad4 = true; quadrant = 4;}
00258                 else return vector<Trajectory>();
00259       }
00260 
00261 
00262         if(staHits.empty()) return vector<Trajectory>();
00263 
00264         if (quad1 && quad2) return vector<Trajectory>();
00265         if (quad1 && quad3) return vector<Trajectory>();
00266         if (quad1 && quad4) return vector<Trajectory>();
00267         if (quad2 && quad3) return vector<Trajectory>();
00268         if (quad2 && quad4) return vector<Trajectory>();
00269         if (quad3 && quad4) return vector<Trajectory>();
00270 
00271         bool doReverse = false;
00272 
00273         if (quad1) doReverse = SlopeSum(staHits); 
00274         if (quad2) doReverse = SlopeSum(staHits); 
00275         if (quad3) doReverse = SlopeSum(staHits); 
00276         if (quad4) doReverse = SlopeSum(staHits); 
00277 
00278 
00279         if(doReverse){
00280         reverse(staHits.begin(),staHits.end());
00281         }
00282 
00283         copy(staHits.begin(),staHits.end(),back_inserter(recHitsForReFit));
00284 
00290 
00291 
00292   if(recHitsForReFit.size() < 2) return vector<Trajectory>();
00293 
00294   bool up = recHitsForReFit.back()->globalPosition().y()>0 ? true : false;
00295   LogTrace(metname) << "Up ? " << up;
00296 
00297   float sumdy = SumDy(recHitsForReFit);
00298 
00299 
00300   PropagationDirection propagationDirection = doReverse ? oppositeToMomentum : alongMomentum;
00301   TrajectoryStateOnSurface firstTSOS = doReverse ? track.outermostMeasurementState() : track.innermostMeasurementState();
00302   unsigned int innerId = doReverse ? track.track().outerDetId() : track.track().innerDetId();
00303 
00304   LogTrace(metname) << "Prop Dir: " << propagationDirection << " FirstId " << innerId << " firstTSOS " << firstTSOS;
00305 
00306   TrajectorySeed seed(PTrajectoryStateOnDet(),TrajectorySeed::recHitContainer(),propagationDirection);
00307 
00308 
00309   if(recHitsForReFit.front()->geographicalId() != DetId(innerId)){
00310     LogTrace(metname)<<"Propagation occurring"<<endl;
00311     firstTSOS = propagator(up, quadrant, sumdy)->propagate(firstTSOS, recHitsForReFit.front()->det()->surface());
00312     LogTrace(metname)<<"Final destination: " << recHitsForReFit.front()->det()->surface().position() << endl;
00313     if(!firstTSOS.isValid()){
00314                 std::cout<<"Propagation error! Dumping RecHits..."<<std::endl;
00315                 
00316                 for(TransientTrackingRecHit::ConstRecHitContainer::const_iterator hit = recHitsForReFit.begin();
00317                 hit !=recHitsForReFit.end(); ++hit){
00318 
00319                 DetId hitId = (*hit)->geographicalId();
00320                 GlobalPoint glbpoint = trackingGeometry()->idToDet(hitId)->position();
00321                 if(hitId.subdetId() == MuonSubdetId::DT)
00322                         std::cout<< glbpoint << " I am DT " << DTWireId(hitId)<<std::endl;
00323                 else if (hitId.subdetId() == MuonSubdetId::CSC )
00324                         std::cout<< glbpoint << " I am CSC " << CSCDetId(hitId)<<std::endl;
00325                 } 
00326 
00327 
00328       LogTrace(metname)<<"Propagation error!"<<endl;
00329       return vector<Trajectory>();
00330     }
00331   }
00332   
00333 
00334   vector<Trajectory> trajectories = fitter(up, quadrant, sumdy)->fit(seed,recHitsForReFit,firstTSOS);
00335   
00336   if(trajectories.empty()){
00337     LogTrace(metname)<<"No Track refitted!"<<endl;
00338     return vector<Trajectory>();
00339   }
00340   
00341   Trajectory trajectoryBW = trajectories.front();
00342     
00343   vector<Trajectory> trajectoriesSM = smoother(up, quadrant, sumdy)->trajectories(trajectoryBW);
00344 
00345   if(trajectoriesSM.empty()){
00346     LogTrace(metname)<<"No Track smoothed!"<<endl;
00347     return vector<Trajectory>();
00348   }
00349   
00350   return trajectoriesSM;
00351 
00352 }
00353 
00354 
00355 
00357 bool TrackTransformerForCosmicMuons::SlopeSum(TransientTrackingRecHit::ConstRecHitContainer tkHits) const{
00358 
00359         bool retval = false;
00360         float y1 = 0 ;
00361         float z1 = 0 ;
00362 
00363         bool first = true;
00364 
00365         std::vector<float> py;
00366         std::vector<float> pz;
00367         int quadrant = -1;
00368         
00369         float sumdy = 0;
00370         float sumdz = 0;
00371 
00372         for(TransientTrackingRecHit::ConstRecHitContainer::const_iterator hit = tkHits.begin();
00373       hit !=tkHits.end(); ++hit){
00374 
00375         DetId hitId = (*hit)->geographicalId();
00376         GlobalPoint glbpoint = trackingGeometry()->idToDet(hitId)->position();
00377         if ( hitId.det() != DetId::Muon || hitId.subdetId() == 3 )continue;
00378 
00379                 float y2 = glbpoint.y();
00380                 float z2 = glbpoint.z();
00381                 float dslope = 0;
00382                 float dy = y2 - y1;
00383                 float dz = z2 - z1;
00384 
00385                 if (y2 > 0 && z2 > 0) quadrant = 1;
00386                 else if (y2 > 0 && z2 < 0) quadrant = 2;
00387                 else if (y2 < 0 && z2 < 0) quadrant = 3;
00388                 else if (y2 < 0 && z2 > 0) quadrant = 4;
00389 
00390 
00391                 if (fabs(dz) > 1e-3) dslope = dy / dz;
00392                 if ( !first) {
00393                         retval+=dslope; 
00394                         sumdy +=dy;
00395                         sumdz +=dz;
00396                 }
00397                 first = false;
00398                 py.push_back(y1);
00399                 pz.push_back(z1);
00400                 y1 = y2;
00401                 z1 = z2;
00402         }
00403 //      std::cout<<"quadrant "<<quadrant;
00404 //      std::cout<<"\tsum dy = "<<sumdy;
00405 //      std::cout<<"\tsum dz = "<<sumdz;
00406 //      std::cout<<std::endl;
00407 
00408 
00409         if ( sumdz < 0) retval = true;
00410         
00411         return retval;
00412 
00413 }
00414 
00415 
00417 float TrackTransformerForCosmicMuons::SumDy(TransientTrackingRecHit::ConstRecHitContainer tkHits) const{
00418 
00419         bool retval = false;
00420         float y1 = 0 ;
00421         float z1 = 0 ;
00422 
00423         bool first = true;
00424 
00425         std::vector<float> py;
00426         std::vector<float> pz;
00427         int quadrant = -1;
00428         
00429         float sumdy = 0;
00430         float sumdz = 0;
00431 
00432 
00433 
00434         for(TransientTrackingRecHit::ConstRecHitContainer::const_iterator hit = tkHits.begin();
00435       hit !=tkHits.end(); ++hit){
00436 
00437         DetId hitId = (*hit)->geographicalId();
00438         GlobalPoint glbpoint = trackingGeometry()->idToDet(hitId)->position();
00439         if ( hitId.det() != DetId::Muon || hitId.subdetId() == 3 )continue;
00440 
00441                 float y2 = glbpoint.y();
00442                 float z2 = glbpoint.z();
00443                 float dslope = 0;
00444                 float dy = y2 - y1;
00445                 float dz = z2 - z1;
00446 
00447                 if (y2 > 0 && z2 > 0) quadrant = 1;
00448                 else if (y2 > 0 && z2 < 0) quadrant = 2;
00449                 else if (y2 < 0 && z2 < 0) quadrant = 3;
00450                 else if (y2 < 0 && z2 > 0) quadrant = 4;
00451 
00452 
00453                 if (fabs(dz) > 1e-3) dslope = dy / dz;
00454                 if ( !first) {
00455                         retval+=dslope; 
00456                         sumdy +=dy;
00457                         sumdz +=dz;
00458                 }
00459                 first = false;
00460                 py.push_back(y1);
00461                 pz.push_back(z1);
00462                 y1 = y2;
00463                 z1 = z2;
00464         }
00465 //      std::cout<<"quadrant "<<quadrant;
00466 //      std::cout<<"\tsum dy = "<<sumdy;
00467 //      std::cout<<"\tsum dz = "<<sumdz;
00468 //      std::cout<<std::endl;
00469 
00470         return sumdy;
00471 }
00472