Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:39 

0001 #include "TrackingTools/TrackRefitter/interface/TrackTransformerForCosmicMuons.h"
0002 
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 
0007 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0008 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0009 
0010 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0011 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0012 #include "TrackingTools/PatternTools/interface/TrajectorySmoother.h"
0013 
0014 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0015 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0016 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0017 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0018 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0019 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0020 
0021 #include "DataFormats/TrackReco/interface/Track.h"
0022 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0023 #include "DataFormats/DetId/interface/DetId.h"
0024 
0025 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0026 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0027 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0028 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0029 
0030 using namespace std;
0031 using namespace edm;
0032 
0033 /// Constructor
0034 TrackTransformerForCosmicMuons::TrackTransformerForCosmicMuons(const ParameterSet& parameterSet,
0035                                                                edm::ConsumesCollector iC)
0036     : theIOpropToken(iC.esConsumes(edm::ESInputTag("", "SmartPropagatorRK"))),
0037       theOIpropToken(iC.esConsumes(edm::ESInputTag("", "SmartPropagatorRKOpposite"))),
0038       thGlobTrackGeoToken(iC.esConsumes()),
0039       theMFToken(iC.esConsumes()),
0040       theIOFitterToken(iC.esConsumes(edm::ESInputTag("", "KFFitterForRefitInsideOut"))),
0041       theOIFitterToken(iC.esConsumes(edm::ESInputTag("", "KFSmootherForRefitInsideOut"))),
0042       theIOSmootherToken(iC.esConsumes(edm::ESInputTag("", "KFFitterForRefitOutsideIn"))),
0043       theOISmootherToken(iC.esConsumes(edm::ESInputTag("", "KFSmootherForRefitOutsideIn"))),
0044       theTkRecHitBuildToken(
0045           iC.esConsumes(edm::ESInputTag("", parameterSet.getParameter<string>("TrackerRecHitBuilder")))),
0046       theMuonRecHitBuildToken(
0047           iC.esConsumes(edm::ESInputTag("", parameterSet.getParameter<string>("MuonRecHitBuilder")))) {
0048   theRPCInTheFit = parameterSet.getParameter<bool>("RefitRPCHits");
0049   theCacheId_TC = theCacheId_GTG = theCacheId_MG = theCacheId_TRH = 0;
0050 }
0051 
0052 /// Destructor
0053 TrackTransformerForCosmicMuons::~TrackTransformerForCosmicMuons() {}
0054 
0055 void TrackTransformerForCosmicMuons::setServices(const EventSetup& setup) {
0056   const std::string metname = "Reco|TrackingTools|TrackTransformer";
0057 
0058   theFitterIO = setup.getHandle(theIOFitterToken);
0059   theFitterOI = setup.getHandle(theOIFitterToken);
0060   theSmootherIO = setup.getHandle(theIOSmootherToken);
0061   theSmootherOI = setup.getHandle(theOISmootherToken);
0062 
0063   unsigned long long newCacheId_TC = setup.get<TrackingComponentsRecord>().cacheIdentifier();
0064 
0065   if (newCacheId_TC != theCacheId_TC) {
0066     LogTrace(metname) << "Tracking Component changed!";
0067     theCacheId_TC = newCacheId_TC;
0068     thePropagatorIO = setup.getHandle(theIOpropToken);
0069     thePropagatorOI = setup.getHandle(theOIpropToken);
0070   }
0071 
0072   // Global Tracking Geometry
0073   unsigned long long newCacheId_GTG = setup.get<GlobalTrackingGeometryRecord>().cacheIdentifier();
0074   if (newCacheId_GTG != theCacheId_GTG) {
0075     LogTrace(metname) << "GlobalTrackingGeometry changed!";
0076     theCacheId_GTG = newCacheId_GTG;
0077     theTrackingGeometry = setup.getHandle(thGlobTrackGeoToken);
0078   }
0079 
0080   // Magfield Field
0081   unsigned long long newCacheId_MG = setup.get<IdealMagneticFieldRecord>().cacheIdentifier();
0082   if (newCacheId_MG != theCacheId_MG) {
0083     LogTrace(metname) << "Magnetic Field changed!";
0084     theCacheId_MG = newCacheId_MG;
0085     theMGField = setup.getHandle(theMFToken);
0086   }
0087 
0088   // Transient Rechit Builders
0089   unsigned long long newCacheId_TRH = setup.get<TransientRecHitRecord>().cacheIdentifier();
0090   if (newCacheId_TRH != theCacheId_TRH) {
0091     theCacheId_TRH = newCacheId_TRH;
0092     LogTrace(metname) << "TransientRecHitRecord changed!";
0093     theTrackerRecHitBuilder = setup.getHandle(theTkRecHitBuildToken);
0094     theMuonRecHitBuilder = setup.getHandle(theMuonRecHitBuildToken);
0095   }
0096 }
0097 
0098 TransientTrackingRecHit::ConstRecHitContainer TrackTransformerForCosmicMuons::getTransientRecHits(
0099     const reco::TransientTrack& track) const {
0100   TransientTrackingRecHit::ConstRecHitContainer tkHits;
0101   TransientTrackingRecHit::ConstRecHitContainer staHits;
0102 
0103   bool printout = false;
0104 
0105   bool quad1 = false;
0106   bool quad2 = false;
0107   bool quad3 = false;
0108   bool quad4 = false;
0109 
0110   for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit)
0111     if ((*hit)->isValid())
0112       if ((*hit)->geographicalId().det() == DetId::Muon) {
0113         if ((*hit)->geographicalId().subdetId() == 3 && !theRPCInTheFit) {
0114           LogTrace("Reco|TrackingTools|TrackTransformer") << "RPC Rec Hit discarged";
0115           continue;
0116         }
0117         staHits.push_back(theMuonRecHitBuilder->build(&**hit));
0118         DetId hitId = staHits.back()->geographicalId();
0119         GlobalPoint glbpoint = trackingGeometry()->idToDet(hitId)->position();
0120         float gpy = glbpoint.y();
0121         float gpz = glbpoint.z();
0122         //      if (gpy != 0 && gpz !=0) slopeSum += gpy / gpz;
0123 
0124         if (gpy > 0 && gpz > 0)
0125           quad1 = true;
0126         else if (gpy > 0 && gpz < 0)
0127           quad2 = true;
0128         else if (gpy < 0 && gpz < 0)
0129           quad3 = true;
0130         else if (gpy < 0 && gpz > 0)
0131           quad4 = true;
0132         else
0133           return tkHits;
0134       }
0135 
0136   if (staHits.empty())
0137     return staHits;
0138 
0139   if (quad1 && quad2)
0140     return tkHits;
0141   if (quad1 && quad3)
0142     return tkHits;
0143   if (quad1 && quad4)
0144     return tkHits;
0145   if (quad2 && quad3)
0146     return tkHits;
0147   if (quad2 && quad4)
0148     return tkHits;
0149   if (quad3 && quad4)
0150     return tkHits;
0151 
0152   bool doReverse = staHits.front()->globalPosition().y() > 0 ? true : false;
0153 
0154   if (quad1)
0155     doReverse = SlopeSum(staHits);
0156   if (quad2)
0157     doReverse = SlopeSum(staHits);
0158   if (quad3)
0159     doReverse = SlopeSum(staHits);
0160   if (quad4)
0161     doReverse = SlopeSum(staHits);
0162   if (doReverse) {
0163     reverse(staHits.begin(), staHits.end());
0164   }
0165 
0166   copy(staHits.begin(), staHits.end(), back_inserter(tkHits));
0167 
0168   ///  if ( quad1 && slopeSum < 0) printout = true;
0169   //  if ( quad2 && slopeSum > 0) printout = true;
0170   ///  if ( quad3 && slopeSum > 0) printout = true;
0171   ///  if ( quad4 && slopeSum < 0) printout = true;
0172   //  swap = printout;
0173 
0174   printout = quad1;
0175 
0176   if (printout)
0177     for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator hit = tkHits.begin(); hit != tkHits.end();
0178          ++hit) {
0179       DetId hitId = (*hit)->geographicalId();
0180       GlobalPoint glbpoint = trackingGeometry()->idToDet(hitId)->position();
0181 
0182       if (hitId.det() == DetId::Muon) {
0183         if (hitId.subdetId() == MuonSubdetId::DT)
0184           LogTrace("TrackFitters") << glbpoint << " I am DT " << DTWireId(hitId);
0185         //      std::cout<< glbpoint << " I am DT " << DTWireId(hitId)<<std::endl;
0186         else if (hitId.subdetId() == MuonSubdetId::CSC)
0187           LogTrace("TrackFitters") << glbpoint << " I am CSC " << CSCDetId(hitId);
0188         //      std::cout<< glbpoint << " I am CSC " << CSCDetId(hitId)<<std::endl;
0189         else if (hitId.subdetId() == MuonSubdetId::RPC)
0190           LogTrace("TrackFitters") << glbpoint << " I am RPC " << RPCDetId(hitId);
0191         else
0192           LogTrace("TrackFitters") << " UNKNOWN MUON HIT TYPE ";
0193       }
0194     }
0195   return tkHits;
0196 }
0197 
0198 /// the refitter used to refit the reco::Track
0199 ESHandle<TrajectoryFitter> TrackTransformerForCosmicMuons::fitter(bool up, int quad, float sumy) const {
0200   if (quad == 1) {
0201     if (sumy < 0)
0202       return theFitterOI;
0203     else
0204       return theFitterIO;
0205   }
0206   if (quad == 2) {
0207     if (sumy < 0)
0208       return theFitterOI;
0209     else
0210       return theFitterIO;
0211   }
0212   if (quad == 3) {
0213     if (sumy > 0)
0214       return theFitterOI;
0215     else
0216       return theFitterIO;
0217   }
0218   if (quad == 4) {
0219     if (sumy > 0)
0220       return theFitterOI;
0221     else
0222       return theFitterIO;
0223   }
0224 
0225   if (up)
0226     return theFitterOI;
0227   else
0228     return theFitterIO;
0229 }
0230 
0231 /// the smoother used to smooth the trajectory which came from the refitting step
0232 ESHandle<TrajectorySmoother> TrackTransformerForCosmicMuons::smoother(bool up, int quad, float sumy) const {
0233   if (quad == 1) {
0234     if (sumy < 0)
0235       return theSmootherOI;
0236     else
0237       return theSmootherIO;
0238   }
0239   if (quad == 2) {
0240     if (sumy < 0)
0241       return theSmootherOI;
0242     else
0243       return theSmootherIO;
0244   }
0245   if (quad == 3) {
0246     if (sumy > 0)
0247       return theSmootherOI;
0248     else
0249       return theSmootherIO;
0250   }
0251   if (quad == 4) {
0252     if (sumy > 0)
0253       return theSmootherOI;
0254     else
0255       return theSmootherIO;
0256   }
0257   if (up)
0258     return theSmootherOI;
0259   else
0260     return theSmootherIO;
0261 }
0262 
0263 ESHandle<Propagator> TrackTransformerForCosmicMuons::propagator(bool up, int quad, float sumy) const {
0264   if (quad == 1) {
0265     if (sumy > 0)
0266       return thePropagatorOI;
0267     else
0268       return thePropagatorIO;
0269   }
0270   if (quad == 2) {
0271     if (sumy > 0)
0272       return thePropagatorOI;
0273     else
0274       return thePropagatorIO;
0275   }
0276   if (quad == 3) {
0277     if (sumy < 0)
0278       return thePropagatorOI;
0279     else
0280       return thePropagatorIO;
0281   }
0282   if (quad == 4) {
0283     if (sumy < 0)
0284       return thePropagatorOI;
0285     else
0286       return thePropagatorIO;
0287   }
0288   if (up)
0289     return thePropagatorIO;
0290   else
0291     return thePropagatorOI;
0292 }
0293 
0294 /// Convert Tracks into Trajectories
0295 vector<Trajectory> TrackTransformerForCosmicMuons::transform(const reco::Track& tr) const {
0296   const std::string metname = "Reco|TrackingTools|TrackTransformer";
0297 
0298   reco::TransientTrack track(tr, magneticField(), trackingGeometry());
0299 
0300   // Build the transient Rechits
0301   TransientTrackingRecHit::ConstRecHitContainer recHitsForReFit;  // = getTransientRecHits(track);
0302   TransientTrackingRecHit::ConstRecHitContainer staHits;
0303 
0304   bool quad1 = false;
0305   bool quad2 = false;
0306   bool quad3 = false;
0307   bool quad4 = false;
0308   int quadrant = 0;
0309 
0310   for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit)
0311     if ((*hit)->isValid())
0312       if ((*hit)->geographicalId().det() == DetId::Muon) {
0313         if ((*hit)->geographicalId().subdetId() == 3 && !theRPCInTheFit) {
0314           LogTrace("Reco|TrackingTools|TrackTransformer") << "RPC Rec Hit discarged";
0315           continue;
0316         }
0317         staHits.push_back(theMuonRecHitBuilder->build(&**hit));
0318         DetId hitId = staHits.back()->geographicalId();
0319         GlobalPoint glbpoint = trackingGeometry()->idToDet(hitId)->position();
0320         float gpy = glbpoint.y();
0321         float gpz = glbpoint.z();
0322         //      if (gpy != 0 && gpz !=0) slopeSum += gpy / gpz;
0323 
0324         if (gpy > 0 && gpz > 0) {
0325           quad1 = true;
0326           quadrant = 1;
0327         } else if (gpy > 0 && gpz < 0) {
0328           quad2 = true;
0329           quadrant = 2;
0330         } else if (gpy < 0 && gpz < 0) {
0331           quad3 = true;
0332           quadrant = 3;
0333         } else if (gpy < 0 && gpz > 0) {
0334           quad4 = true;
0335           quadrant = 4;
0336         } else
0337           return vector<Trajectory>();
0338       }
0339 
0340   if (staHits.empty())
0341     return vector<Trajectory>();
0342 
0343   if (quad1 && quad2)
0344     return vector<Trajectory>();
0345   if (quad1 && quad3)
0346     return vector<Trajectory>();
0347   if (quad1 && quad4)
0348     return vector<Trajectory>();
0349   if (quad2 && quad3)
0350     return vector<Trajectory>();
0351   if (quad2 && quad4)
0352     return vector<Trajectory>();
0353   if (quad3 && quad4)
0354     return vector<Trajectory>();
0355 
0356   bool doReverse = false;
0357 
0358   if (quad1)
0359     doReverse = SlopeSum(staHits);
0360   if (quad2)
0361     doReverse = SlopeSum(staHits);
0362   if (quad3)
0363     doReverse = SlopeSum(staHits);
0364   if (quad4)
0365     doReverse = SlopeSum(staHits);
0366 
0367   if (doReverse) {
0368     reverse(staHits.begin(), staHits.end());
0369   }
0370 
0371   copy(staHits.begin(), staHits.end(), back_inserter(recHitsForReFit));
0372 
0373   ///
0374   ///
0375   ///
0376   ///
0377   ///
0378 
0379   if (recHitsForReFit.size() < 2)
0380     return vector<Trajectory>();
0381 
0382   bool up = recHitsForReFit.back()->globalPosition().y() > 0 ? true : false;
0383   LogTrace(metname) << "Up ? " << up;
0384 
0385   float sumdy = SumDy(recHitsForReFit);
0386 
0387   PropagationDirection propagationDirection = doReverse ? oppositeToMomentum : alongMomentum;
0388   TrajectoryStateOnSurface firstTSOS =
0389       doReverse ? track.outermostMeasurementState() : track.innermostMeasurementState();
0390   unsigned int innerId = doReverse ? track.track().outerDetId() : track.track().innerDetId();
0391 
0392   LogTrace(metname) << "Prop Dir: " << propagationDirection << " FirstId " << innerId << " firstTSOS " << firstTSOS;
0393 
0394   TrajectorySeed seed({}, {}, propagationDirection);
0395 
0396   if (recHitsForReFit.front()->geographicalId() != DetId(innerId)) {
0397     LogTrace(metname) << "Propagation occurring" << endl;
0398     firstTSOS = propagator(up, quadrant, sumdy)->propagate(firstTSOS, recHitsForReFit.front()->det()->surface());
0399     LogTrace(metname) << "Final destination: " << recHitsForReFit.front()->det()->surface().position() << endl;
0400     if (!firstTSOS.isValid()) {
0401       std::cout << "Propagation error! Dumping RecHits..." << std::endl;
0402 
0403       for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator hit = recHitsForReFit.begin();
0404            hit != recHitsForReFit.end();
0405            ++hit) {
0406         DetId hitId = (*hit)->geographicalId();
0407         GlobalPoint glbpoint = trackingGeometry()->idToDet(hitId)->position();
0408         if (hitId.subdetId() == MuonSubdetId::DT)
0409           std::cout << glbpoint << " I am DT " << DTWireId(hitId) << std::endl;
0410         else if (hitId.subdetId() == MuonSubdetId::CSC)
0411           std::cout << glbpoint << " I am CSC " << CSCDetId(hitId) << std::endl;
0412       }
0413 
0414       LogTrace(metname) << "Propagation error!" << endl;
0415       return vector<Trajectory>();
0416     }
0417   }
0418 
0419   vector<Trajectory>&& trajectories = fitter(up, quadrant, sumdy)->fit(seed, recHitsForReFit, firstTSOS);
0420 
0421   if (trajectories.empty()) {
0422     LogTrace(metname) << "No Track refitted!" << endl;
0423     return vector<Trajectory>();
0424   }
0425 
0426   Trajectory const& trajectoryBW = trajectories.front();
0427 
0428   vector<Trajectory> trajectoriesSM = smoother(up, quadrant, sumdy)->trajectories(trajectoryBW);
0429 
0430   if (trajectoriesSM.empty()) {
0431     LogTrace(metname) << "No Track smoothed!" << endl;
0432   }
0433 
0434   return trajectoriesSM;
0435 }
0436 
0437 ///decide if the track should be reversed
0438 bool TrackTransformerForCosmicMuons::SlopeSum(const TransientTrackingRecHit::ConstRecHitContainer& tkHits) const {
0439   bool retval = false;
0440   float y1 = 0;
0441   float z1 = 0;
0442 
0443   bool first = true;
0444 
0445   std::vector<float> py;
0446   std::vector<float> pz;
0447   //    int quadrant = -1;
0448 
0449   float sumdz = 0;
0450 
0451   for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator hit = tkHits.begin(); hit != tkHits.end(); ++hit) {
0452     DetId hitId = (*hit)->geographicalId();
0453     GlobalPoint glbpoint = trackingGeometry()->idToDet(hitId)->position();
0454     if (hitId.det() != DetId::Muon || hitId.subdetId() == 3)
0455       continue;
0456 
0457     float y2 = glbpoint.y();
0458     float z2 = glbpoint.z();
0459     float dslope = 0;
0460     float dy = y2 - y1;
0461     float dz = z2 - z1;
0462 
0463     //      if (y2 > 0 && z2 > 0) quadrant = 1;
0464     //      else if (y2 > 0 && z2 < 0) quadrant = 2;
0465     //      else if (y2 < 0 && z2 < 0) quadrant = 3;
0466     //      else if (y2 < 0 && z2 > 0) quadrant = 4;
0467 
0468     if (fabs(dz) > 1e-3)
0469       dslope = dy / dz;
0470     if (!first) {
0471       retval += dslope;
0472       sumdz += dz;
0473     }
0474     first = false;
0475     py.push_back(y1);
0476     pz.push_back(z1);
0477     y1 = y2;
0478     z1 = z2;
0479   }
0480   //    std::cout<<"quadrant "<<quadrant;
0481   //    std::cout<<"\tsum dy = "<<sumdy;
0482   //    std::cout<<"\tsum dz = "<<sumdz;
0483   //    std::cout<<std::endl;
0484 
0485   if (sumdz < 0)
0486     retval = true;
0487 
0488   return retval;
0489 }
0490 
0491 ///decide if the track should be reversed
0492 float TrackTransformerForCosmicMuons::SumDy(const TransientTrackingRecHit::ConstRecHitContainer& tkHits) const {
0493   float y1 = 0;
0494   float z1 = 0;
0495 
0496   bool first = true;
0497 
0498   std::vector<float> py;
0499   std::vector<float> pz;
0500   //    int quadrant = -1;
0501 
0502   float sumdy = 0;
0503 
0504   for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator hit = tkHits.begin(); hit != tkHits.end(); ++hit) {
0505     DetId hitId = (*hit)->geographicalId();
0506     GlobalPoint glbpoint = trackingGeometry()->idToDet(hitId)->position();
0507     if (hitId.det() != DetId::Muon || hitId.subdetId() == 3)
0508       continue;
0509 
0510     float y2 = glbpoint.y();
0511     float z2 = glbpoint.z();
0512     float dy = y2 - y1;
0513 
0514     //      if (y2 > 0 && z2 > 0) quadrant = 1;
0515     //      else if (y2 > 0 && z2 < 0) quadrant = 2;
0516     //      else if (y2 < 0 && z2 < 0) quadrant = 3;
0517     //      else if (y2 < 0 && z2 > 0) quadrant = 4;
0518 
0519     if (!first) {
0520       sumdy += dy;
0521     }
0522     first = false;
0523     py.push_back(y1);
0524     pz.push_back(z1);
0525     y1 = y2;
0526     z1 = z2;
0527   }
0528   //    std::cout<<"quadrant "<<quadrant;
0529   //    std::cout<<"\tsum dy = "<<sumdy;
0530   //    std::cout<<"\tsum dz = "<<sumdz;
0531   //    std::cout<<std::endl;
0532 
0533   return sumdy;
0534 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!