TrackRefitter/src/TrackTransformer.cc

0001 #include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"
0002
0003 #include "FWCore/Framework/interface/ConsumesCollector.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0009 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0010
0011 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0012 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0013 #include "TrackingTools/PatternTools/interface/TrajectorySmoother.h"
0014
0015 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0016 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0017 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0018 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0019 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0020 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0021
0022 #include "DataFormats/TrackReco/interface/Track.h"
0023 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0024 #include "DataFormats/DetId/interface/DetId.h"
0025
0026 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0027
0028 using namespace std;
0029 using namespace edm;
0030
0031 /// Constructor
0032 TrackTransformer::TrackTransformer(const ParameterSet& parameterSet, edm::ConsumesCollector& iC)
0033     : theRPCInTheFit(parameterSet.getParameter<bool>("RefitRPCHits")),
0034       theDoPredictionsOnly(parameterSet.getParameter<bool>("DoPredictionsOnly")),
0035       theRefitDirection(parameterSet.getParameter<string>("RefitDirection")),
0036       theFitterName(parameterSet.getParameter<string>("Fitter")),
0037       theSmootherName(parameterSet.getParameter<string>("Smoother")),
0038       thePropagatorName(parameterSet.getParameter<string>("Propagator")),
0039       theTrackerRecHitBuilderName(parameterSet.getParameter<string>("TrackerRecHitBuilder")),
0040       theMuonRecHitBuilderName(parameterSet.getParameter<string>("MuonRecHitBuilder")),
0041       theMTDRecHitBuilderName(parameterSet.getParameter<string>("MTDRecHitBuilder")) {
0042   theTrackingGeometryToken = iC.esConsumes();
0043   theMGFieldToken = iC.esConsumes();
0044   theFitterToken = iC.esConsumes(edm::ESInputTag("", theFitterName));
0045   theSmootherToken = iC.esConsumes(edm::ESInputTag("", theSmootherName));
0046   thePropagatorToken = iC.esConsumes(edm::ESInputTag("", thePropagatorName));
0047   theTrackerRecHitBuilderToken = iC.esConsumes(edm::ESInputTag("", theTrackerRecHitBuilderName));
0048   theMuonRecHitBuilderToken = iC.esConsumes(edm::ESInputTag("", theMuonRecHitBuilderName));
0049   theMTDRecHitBuilderToken = iC.esConsumes(edm::ESInputTag("", theMTDRecHitBuilderName));
0050 }
0051
0052 /// Destructor
0053 TrackTransformer::~TrackTransformer() {}
0054
0055 void TrackTransformer::fillPSetDescription(edm::ParameterSetDescription& desc,
0056                                            bool doPredictionsOnly,
0057                                            const std::string& fitter,
0058                                            const std::string& smoother,
0059                                            const std::string& propagator,
0060                                            const std::string& refitDirection,
0061                                            bool refitRPCHits,
0062                                            const std::string& trackerRecHitBuilder,
0063                                            const std::string& muonRecHitBuilder,
0064                                            const std::string& mtdRecHitBuilder) {
0065   desc.add<bool>("DoPredictionsOnly", doPredictionsOnly);
0066   desc.add<std::string>("Fitter", fitter);
0067   desc.add<std::string>("Smoother", smoother);
0068   desc.add<std::string>("Propagator", propagator);
0069   desc.add<std::string>("RefitDirection", refitDirection);
0070   desc.add<bool>("RefitRPCHits", refitRPCHits);
0071   desc.add<std::string>("TrackerRecHitBuilder", trackerRecHitBuilder);
0072   desc.add<std::string>("MuonRecHitBuilder", muonRecHitBuilder);
0073   desc.add<std::string>("MTDRecHitBuilder", mtdRecHitBuilder);
0074 }
0075
0076 void TrackTransformer::setServices(const EventSetup& setup) {
0077   const std::string metname = "Reco|TrackingTools|TrackTransformer";
0078
0079   theFitter = setup.getData(theFitterToken).clone();
0080   theSmoother.reset(setup.getData(theSmootherToken).clone());
0081
0082   thePropagator = setup.getHandle(thePropagatorToken);
0083
0084   // Global Tracking Geometry
0085   theTrackingGeometry = setup.getHandle(theTrackingGeometryToken);
0086
0087   // Magfield Field
0088   theMGField = setup.getHandle(theMGFieldToken);
0089
0090   // Transient Rechit Builders
0091   unsigned long long newCacheId_TRH = setup.get<TransientRecHitRecord>().cacheIdentifier();
0092   if (newCacheId_TRH != theCacheId_TRH) {
0093     theCacheId_TRH = newCacheId_TRH;
0094     LogTrace(metname) << "TransientRecHitRecord changed!";
0095     theTrackerRecHitBuilder = &setup.getData(theTrackerRecHitBuilderToken);
0096     theMuonRecHitBuilder = setup.getHandle(theMuonRecHitBuilderToken);
0097     theMTDRecHitBuilder = setup.getHandle(theMTDRecHitBuilderToken);
0098     theMtdAvailable = theMTDRecHitBuilder.isValid();
0099     hitCloner = static_cast<TkTransientTrackingRecHitBuilder const*>(theTrackerRecHitBuilder)->cloner();
0100   }
0101   theFitter->setHitCloner(&hitCloner);
0102   theSmoother->setHitCloner(&hitCloner);
0103 }
0104
0105 vector<Trajectory> TrackTransformer::transform(const reco::TrackRef& track) const { return transform(*track); }
0106
0107 TransientTrackingRecHit::ConstRecHitContainer TrackTransformer::getTransientRecHits(
0108     const reco::TransientTrack& track) const {
0109   TransientTrackingRecHit::ConstRecHitContainer result;
0110   auto tkbuilder = static_cast<TkTransientTrackingRecHitBuilder const*>(theTrackerRecHitBuilder);
0111
0112   for (auto hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
0113     if ((*hit)->isValid()) {
0114       if ((*hit)->geographicalId().det() == DetId::Tracker) {
0115         result.emplace_back((**hit).cloneForFit(*tkbuilder->geometry()->idToDet((**hit).geographicalId())));
0116       } else if ((*hit)->geographicalId().det() == DetId::Muon) {
0117         if ((*hit)->geographicalId().subdetId() == 3 && !theRPCInTheFit) {
0118           LogTrace("Reco|TrackingTools|TrackTransformer") << "RPC Rec Hit discarged";
0119           continue;
0120         }
0121         result.push_back(theMuonRecHitBuilder->build(&**hit));
0122       } else if ((*hit)->geographicalId().det() == DetId::Forward && (*hit)->geographicalId().subdetId() == FastTime) {
0123         if (theMtdAvailable)
0124           result.push_back(theMTDRecHitBuilder->build(&**hit));
0125         else
0126           throw cms::Exception("TrackTransformer") << "MTD hit encountered but MTD not available!";
0127       }
0128     }
0129   }
0130
0131   return result;
0132 }
0133
0134 // FIXME: check this method!
0135 RefitDirection::GeometricalDirection TrackTransformer::checkRecHitsOrdering(
0136     TransientTrackingRecHit::ConstRecHitContainer const& recHits) const {
0137   if (!recHits.empty()) {
0138     GlobalPoint first = trackingGeometry()->idToDet(recHits.front()->geographicalId())->position();
0139     GlobalPoint last = trackingGeometry()->idToDet(recHits.back()->geographicalId())->position();
0140
0141     // maybe perp2?
0142     auto rFirst = first.mag2();
0143     auto rLast = last.mag2();
0144     if (rFirst < rLast)
0145       return RefitDirection::insideOut;
0146     if (rFirst > rLast)
0147       return RefitDirection::outsideIn;
0148   }
0149   LogDebug("Reco|TrackingTools|TrackTransformer") << "Impossible to determine the rechits order" << endl;
0150   return RefitDirection::undetermined;
0151 }
0152
0153 /// Convert Tracks into Trajectories
0154 vector<Trajectory> TrackTransformer::transform(const reco::Track& newTrack) const {
0155   const std::string metname = "Reco|TrackingTools|TrackTransformer";
0156
0157   reco::TransientTrack track(newTrack, magneticField(), trackingGeometry());
0158
0159   auto recHitsForReFit = getTransientRecHits(track);
0160   return transform(track, recHitsForReFit);
0161 }
0162
0163 /// Convert Tracks into Trajectories with a given set of hits
0164 vector<Trajectory> TrackTransformer::transform(const reco::TransientTrack& track,
0165                                                TransientTrackingRecHit::ConstRecHitContainer& recHitsForReFit) const {
0166   const std::string metname = "Reco|TrackingTools|TrackTransformer";
0167
0168   if (recHitsForReFit.size() < 2)
0169     return vector<Trajectory>();
0170
0171   // 8 cases are foreseen:
0172   // [RH = rec hit order, P = momentum dir, FD = fit direction. IO/OI = inside-out/outside-in, AM/OM = along momentum/opposite to momentum]
0173   // (1) RH IO | P IO | FD AM  ---> Start from IN
0174   // (2) RH IO | P IO | FD OM  ---> Reverse RH and start from OUT
0175   // (3) RH IO | P OI | FD AM  ---> Reverse RH and start from IN
0176   // (4) RH IO | P OI | FD OM  ---> Start from OUT
0177   // (5) RH OI | P IO | FD AM  ---> Reverse RH and start from IN
0178   // (6) RH OI | P IO | FD OM  ---> Start from OUT
0179   // (7) RH OI | P OI | FD AM  ---> Start from IN
0180   // (8) RH OI | P OI | FD OM  ---> Reverse RH and start from OUT
0181   //
0182   // *** Rules: ***
0183   // -A- If RH-FD agree (IO-AM,OI-OM) do not reverse the RH
0184   // -B- If FD along momentum start from innermost state, otherwise use outermost
0185
0186   // Other special cases can be handled:
0187   // (1 bis) RH IO | P IO | GFD IO => FD AM  ---> Start from IN
0188   // (2 bis) RH IO | P IO | GFD OI => FD OM  ---> Reverse RH and start from OUT
0189   // (3 bis) RH IO | P OI | GFD OI => FD AM  ---> Reverse RH and start from OUT
0190   // (4 bis) RH IO | P OI | GFD IO => FD OM  ---> Start from IN
0191   // (5 bis) RH OI | P IO | GFD IO => FD AM  ---> Reverse RH and start from IN
0192   // (6 bis) RH OI | P IO | GFD OI => FD OM  ---> Start from OUT
0193   // (7 bis) RH OI | P OI | GFD OI => FD AM  ---> Start from OUT
0194   // (8 bis) RH OI | P OI | GFD IO => FD OM  ---> Reverse RH and start from IN
0195   //
0196   // *** Additional rule: ***
0197   // -A0- If P and GFD agree, then FD is AM otherwise is OM
0198   // -A00- rechit must be ordered as GFD in order to handle the case of cosmics
0199   // -B0- The starting state is decided by GFD
0200
0201   // Determine the RH order
0202   RefitDirection::GeometricalDirection recHitsOrder =
0203       checkRecHitsOrdering(recHitsForReFit);  // FIXME change nome of the *type*  --> RecHit order!
0204   LogTrace(metname) << "RH order (0-insideOut, 1-outsideIn): " << recHitsOrder;
0205
0206   PropagationDirection propagationDirection = theRefitDirection.propagationDirection();
0207
0208   // Apply rule -A0-
0209   if (propagationDirection == anyDirection) {
0210     GlobalVector momentum = track.innermostMeasurementState().globalMomentum();
0211     GlobalVector position = track.innermostMeasurementState().globalPosition() - GlobalPoint(0, 0, 0);
0212     RefitDirection::GeometricalDirection p = (momentum.x() * position.x() > 0 || momentum.y() * position.y() > 0)
0213                                                  ? RefitDirection::insideOut
0214                                                  : RefitDirection::outsideIn;
0215
0216     propagationDirection = p == theRefitDirection.geometricalDirection() ? alongMomentum : oppositeToMomentum;
0217     LogTrace(metname) << "P  (0-insideOut, 1-outsideIn): " << p;
0218     LogTrace(metname) << "FD (0-OM, 1-AM, 2-ANY): " << propagationDirection;
0219   }
0220   // -A0-
0221
0222   // Apply rule -A-
0223   if (theRefitDirection.propagationDirection() != anyDirection) {
0224     if ((recHitsOrder == RefitDirection::insideOut && propagationDirection == oppositeToMomentum) ||
0225         (recHitsOrder == RefitDirection::outsideIn && propagationDirection == alongMomentum))
0226       reverse(recHitsForReFit.begin(), recHitsForReFit.end());
0227   }
0228   // -A-
0229   // Apply rule -A00-
0230   else {
0231     // reorder the rechit as defined in theRefitDirection.geometricalDirection();
0232     if (theRefitDirection.geometricalDirection() != recHitsOrder)
0233       reverse(recHitsForReFit.begin(), recHitsForReFit.end());
0234   }
0235   // -A00-
0236
0237   // Apply rule -B-
0238   TrajectoryStateOnSurface firstTSOS = track.innermostMeasurementState();
0239   unsigned int innerId = track.track().innerDetId();
0240   if (theRefitDirection.propagationDirection() != anyDirection) {
0241     if (propagationDirection == oppositeToMomentum) {
0242       innerId = track.track().outerDetId();
0243       firstTSOS = track.outermostMeasurementState();
0244     }
0245   } else {  // if(theRefitDirection.propagationDirection() == anyDirection)
0246     // Apply rule -B0-
0247     if (theRefitDirection.geometricalDirection() == RefitDirection::outsideIn) {
0248       innerId = track.track().outerDetId();
0249       firstTSOS = track.outermostMeasurementState();
0250     }
0251     // -B0-
0252   }
0253   // -B-
0254
0255   if (!firstTSOS.isValid()) {
0256     LogTrace(metname) << "Error wrong initial state!" << endl;
0257     return vector<Trajectory>();
0258   }
0259
0260   TrajectorySeed seed({}, {}, propagationDirection);
0261
0262   if (recHitsForReFit.front()->geographicalId() != DetId(innerId)) {
0263     LogTrace(metname) << "Propagation occured" << endl;
0264     firstTSOS = propagator()->propagate(firstTSOS, recHitsForReFit.front()->det()->surface());
0265     if (!firstTSOS.isValid()) {
0266       LogTrace(metname) << "Propagation error!" << endl;
0267       return vector<Trajectory>();
0268     }
0269   }
0270
0271   if (theDoPredictionsOnly) {
0272     Trajectory aTraj(seed, propagationDirection);
0273     TrajectoryStateOnSurface predTSOS = firstTSOS;
0274     for (auto const& hit : recHitsForReFit) {
0275       predTSOS = propagator()->propagate(predTSOS, hit->det()->surface());
0276       if (predTSOS.isValid())
0277         aTraj.push(TrajectoryMeasurement(predTSOS, hit));
0278     }
0279     return vector<Trajectory>(1, aTraj);
0280   }
0281
0282   auto const& trajectories = theFitter->fit(seed, recHitsForReFit, firstTSOS);
0283
0284   if (trajectories.empty()) {
0285     LogTrace(metname) << "No Track refitted!" << endl;
0286     return trajectories;
0287   }
0288
0289   auto const& trajectoryBW = trajectories.front();
0290
0291   auto const& trajectoriesSM = theSmoother->trajectories(trajectoryBW);
0292
0293   if (trajectoriesSM.empty()) {
0294     LogTrace(metname) << "No Track smoothed!" << endl;
0295   }
0296
0297   return trajectoriesSM;
0298 }