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
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
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
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
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
00156
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
00171 else if (hitId.subdetId() == MuonSubdetId::CSC )
00172 LogTrace("TrackFitters") << glbpoint << " I am CSC " << CSCDetId(hitId);
00173
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
00224 TransientTrackingRecHit::ConstRecHitContainer recHitsForReFit;
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
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
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
00380
00381
00382
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
00398
00399
00400
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
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
00442
00443
00444
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
00460
00461
00462
00463
00464 return sumdy;
00465 }
00466