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
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
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
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
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
00162
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
00177 else if (hitId.subdetId() == MuonSubdetId::CSC )
00178 LogTrace("TrackFitters") << glbpoint << " I am CSC " << CSCDetId(hitId);
00179
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
00230 TransientTrackingRecHit::ConstRecHitContainer recHitsForReFit;
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
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
00404
00405
00406
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
00466
00467
00468
00469
00470 return sumdy;
00471 }
00472