CMS 3D CMS Logo

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