CkfPattern/src/BaseCkfTrajectoryBuilder.cc

0001 #include "RecoTracker/CkfPattern/interface/BaseCkfTrajectoryBuilder.h"
0002
0003 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0004 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0005 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0006 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
0007 #include "TrackingTools/TrajectoryFiltering/interface/TrajectoryFilter.h"
0008
0009 #include "TrackingTools/DetLayers/interface/NavigationSchool.h"
0010 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0011 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0012 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0013 #include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h"
0014 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
0015 #include "TrackingTools/TrajectoryFiltering/interface/TrajectoryFilterFactory.h"
0016
0017 #include "TrackingTools/PatternTools/interface/TempTrajectory.h"
0018 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0019 #include "TrackingTools/TrajectoryFiltering/interface/TrajectoryFilterFactory.h"
0020
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0024 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0025
0026 BaseCkfTrajectoryBuilder::BaseCkfTrajectoryBuilder(const edm::ParameterSet& conf,
0027                                                    edm::ConsumesCollector iC,
0028                                                    std::unique_ptr<TrajectoryFilter> filter,
0029                                                    std::unique_ptr<TrajectoryFilter> inOutFilter)
0030     : theSeedAs5DHit(conf.getParameter<bool>("seedAs5DHit")),
0031       theFilter(std::move(filter)),
0032       theInOutFilter(std::move(inOutFilter)),
0033       theUpdatorToken(iC.esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("updator")))),
0034       thePropagatorAlongToken(iC.esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("propagatorAlong")))),
0035       thePropagatorOppositeToken(
0036           iC.esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("propagatorOpposite")))),
0037       theEstimatorToken(iC.esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("estimator")))),
0038       theRecHitBuilderToken(iC.esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("TTRHBuilder")))) {
0039   if (conf.exists("clustersToSkip"))
0040     edm::LogError("BaseCkfTrajectoryBuilder")
0041         << "ERROR: " << typeid(*this).name() << " has a clustersToSkip parameter set";
0042 }
0043
0044 BaseCkfTrajectoryBuilder::~BaseCkfTrajectoryBuilder() {}
0045
0046 void BaseCkfTrajectoryBuilder::fillPSetDescription(edm::ParameterSetDescription& iDesc) {
0047   iDesc.add<bool>("seedAs5DHit", false);
0048   iDesc.add<std::string>("updator", "KFUpdator");
0049   iDesc.add<std::string>("propagatorAlong", "PropagatorWithMaterial");
0050   iDesc.add<std::string>("propagatorOpposite", "PropagatorWithMaterialOpposite");
0051   iDesc.add<std::string>("estimator", "Chi2");
0052   iDesc.add<std::string>("TTRHBuilder", "WithTrackAngle");
0053 }
0054
0055 std::unique_ptr<TrajectoryFilter> BaseCkfTrajectoryBuilder::createTrajectoryFilter(const edm::ParameterSet& pset,
0056                                                                                    edm::ConsumesCollector& iC) {
0057   return TrajectoryFilterFactory::get()->create(pset.getParameter<std::string>("ComponentType"), pset, iC);
0058 }
0059
0060 #include "RecoTracker/TransientTrackingRecHit/interface/TRecHit5DParamConstraint.h"
0061 void BaseCkfTrajectoryBuilder::seedMeasurements(const TrajectorySeed& seed, TempTrajectory& result, bool as5D) const {
0062   PTrajectoryStateOnDet pState(seed.startingState());
0063   const GeomDet* gdet = theMeasurementTracker->geomTracker()->idToDet(pState.detId());
0064   TSOS outerState =
0065       trajectoryStateTransform::transientState(pState, &(gdet->surface()), forwardPropagator(seed)->magneticField());
0066
0067   if (as5D) {
0068     TrackingRecHit::RecHitPointer recHit(new TRecHit5DParamConstraint(*gdet, outerState));
0069     TSOS invalidState(gdet->surface());
0070     auto hitLayer = theMeasurementTracker->geometricSearchTracker()->detLayer(pState.detId());
0071     result.emplace(invalidState, outerState, recHit, 0, hitLayer);
0072     return;
0073   }
0074
0075   for (auto ihit = seed.recHits().begin(); ihit != seed.recHits().end(); ihit++) {
0076     TrackingRecHit::RecHitPointer recHit = ihit->cloneSH();
0077     const GeomDet* hitGeomDet = recHit->det();
0078
0079     const DetLayer* hitLayer = theMeasurementTracker->geometricSearchTracker()->detLayer(ihit->geographicalId());
0080
0081     TSOS invalidState(hitGeomDet->surface());
0082     if (ihit == seed.recHits().end() - 1) {
0083       // the seed trajectory state should correspond to this hit
0084       if (&gdet->surface() != &hitGeomDet->surface()) {
0085         edm::LogError("CkfPattern")
0086             << "CkfTrajectoryBuilder error: the seed state is not on the surface of the detector of the last seed hit";
0087         return;  // FIXME: should throw exception
0088       }
0089
0090       //TSOS updatedState = outerstate;
0091       result.emplace(invalidState, outerState, recHit, 0, hitLayer);
0092     } else {
0093       TSOS innerState = backwardPropagator(seed)->propagate(outerState, hitGeomDet->surface());
0094
0095       // try to recover if propagation failed
0096       if UNLIKELY (!innerState.isValid())
0097         innerState = trajectoryStateTransform::transientState(
0098             pState, &(hitGeomDet->surface()), forwardPropagator(seed)->magneticField());
0099
0100       if (innerState.isValid()) {
0101         TSOS innerUpdated = theUpdator->update(innerState, *recHit);
0102         result.emplace(invalidState, innerUpdated, recHit, 0, hitLayer);
0103       }
0104     }
0105   }
0106
0107   // method for debugging
0108   // fix somehow
0109   // fillSeedHistoDebugger(result.begin(),result.end());
0110 }
0111
0112 TempTrajectory BaseCkfTrajectoryBuilder::createStartingTrajectory(const TrajectorySeed& seed) const {
0113   TempTrajectory result(seed.direction(), seed.nHits());
0114   seedMeasurements(seed, result, theSeedAs5DHit);
0115
0116   LogDebug("CkfPattern") << " initial trajectory from the seed: " << PrintoutHelper::dumpCandidate(result, true);
0117
0118   return result;
0119 }
0120
0121 bool BaseCkfTrajectoryBuilder::toBeContinued(TempTrajectory& traj, bool inOut) const {
0122   if UNLIKELY (traj.measurements().size() > 400) {
0123     edm::LogError("BaseCkfTrajectoryBuilder_InfiniteLoop");
0124     LogTrace("BaseCkfTrajectoryBuilder_InfiniteLoop")
0125         << "Cropping Track After 400 Measurements:\n"
0126         << "   Last predicted state: " << traj.lastMeasurement().predictedState() << "\n"
0127         << "   Last layer subdetector: " << (traj.lastLayer() ? traj.lastLayer()->subDetector() : -1) << "\n"
0128         << "   Found hits: " << traj.foundHits() << ", lost hits: " << traj.lostHits() << "\n\n";
0129     return false;
0130   }
0131   // Called after each new hit is added to the trajectory, to see if it is
0132   // worth continuing to build this track candidate.
0133   //
0134   // When a sufficient amount of measurements are made,
0135   // ensure that an infinite loop is not created (CMSHLT-3557).
0136   // Avoid hit-pair structures as last = last-2, and last-1 = last-3,
0137   // where last refers to measurements.
0138   //
0139   const TempTrajectory::DataContainer tms = traj.measurements();
0140   TempTrajectory::DataContainer::const_iterator tm = tms.begin();
0141
0142   // Ensure at sufficient amount of measurements before checking for loops
0143   if (traj.measurements().size() > 15) {
0144     TrackingRecHit::RecHitPointer lastHit = tm->recHit();
0145     ++tm;
0146     TrackingRecHit::RecHitPointer last2Hit = tm->recHit();
0147     ++tm;
0148     TrackingRecHit::RecHitPointer last3Hit = tm->recHit();
0149     ++tm;
0150     TrackingRecHit::RecHitPointer last4Hit = tm->recHit();
0151     if (lastHit->geographicalId() == last3Hit->geographicalId() &&
0152         last2Hit->geographicalId() == last4Hit->geographicalId() &&
0153         (lastHit->geographicalId().rawId() == 0 || last2Hit->geographicalId().rawId() == 0)) {
0154       LogDebug("CkfPattern") << "Loop pattern found in last recHits\n" << PrintoutHelper::dumpMeasurements(tms);
0155
0156       return false;
0157     }
0158   }
0159   if (inOut) {
0160     // if (theInOutFilter == 0) edm::LogError("CkfPattern") << "CkfTrajectoryBuilder error: trying to use dedicated filter for in-out tracking phase, when none specified";
0161     return theInOutFilter->toBeContinued(traj);
0162   } else {
0163     return theFilter->toBeContinued(traj);
0164   }
0165 }
0166
0167 bool BaseCkfTrajectoryBuilder::qualityFilter(const TempTrajectory& traj, bool inOut) const {
0168   // Called after building a trajectory is completed, to see if it is good enough
0169   // to keep.
0170   if (inOut) {
0171     // if (theInOutFilter == 0) edm::LogError("CkfPattern") << "CkfTrajectoryBuilder error: trying to use dedicated filter for in-out tracking phase, when none specified";
0172     return theInOutFilter->qualityFilter(traj);
0173   } else {
0174     return theFilter->qualityFilter(traj);
0175   }
0176 }
0177
0178 void BaseCkfTrajectoryBuilder::addToResult(std::shared_ptr<const TrajectorySeed> const& seed,
0179                                            TempTrajectory& tmptraj,
0180                                            TrajectoryContainer& result,
0181                                            bool inOut) const {
0182   // quality check
0183   if (!qualityFilter(tmptraj, inOut))
0184     return;
0185   Trajectory traj = tmptraj.toTrajectory();
0186   traj.setSharedSeed(seed);
0187   // discard latest dummy measurements
0188   while (!traj.empty() && !traj.lastMeasurement().recHit()->isValid())
0189     traj.pop();
0190   LogDebug("CkfPattern") << inOut << "=inOut option. pushing a Trajectory with: " << traj.foundHits() << " found hits. "
0191                          << traj.lostHits()
0192                          << " lost hits. Popped :" << (tmptraj.measurements().size()) - (traj.measurements().size())
0193                          << " hits.";
0194   result.push_back(std::move(traj));
0195 }
0196
0197 void BaseCkfTrajectoryBuilder::addToResult(TempTrajectory const& tmptraj,
0198                                            TempTrajectoryContainer& result,
0199                                            bool inOut) const {
0200   // quality check
0201   if (!qualityFilter(tmptraj, inOut))
0202     return;
0203   // discard latest dummy measurements
0204   TempTrajectory traj = tmptraj;
0205   while (!traj.empty() && !traj.lastMeasurement().recHit()->isValid())
0206     traj.pop();
0207   LogDebug("CkfPattern") << inOut << "=inOut option. pushing a TempTrajectory with: " << traj.foundHits()
0208                          << " found hits. " << traj.lostHits()
0209                          << " lost hits. Popped :" << (tmptraj.measurements().size()) - (traj.measurements().size())
0210                          << " hits.";
0211   result.push_back(std::move(traj));
0212 }
0213
0214 void BaseCkfTrajectoryBuilder::moveToResult(TempTrajectory&& traj, TempTrajectoryContainer& result, bool inOut) const {
0215   // quality check
0216   if (!qualityFilter(traj, inOut))
0217     return;
0218   // discard latest dummy measurements
0219   while (!traj.empty() && !traj.lastMeasurement().recHitR().isValid())
0220     traj.pop();
0221   LogDebug("CkfPattern") << inOut << "=inOut option. pushing a TempTrajectory with: " << traj.foundHits()
0222                          << " found hits. " << traj.lostHits();
0223   //            <<" lost hits. Popped :"<<(ttraj.measurements().size())-(traj.measurements().size())<<" hits.";
0224   result.push_back(std::move(traj));
0225 }
0226
0227 BaseCkfTrajectoryBuilder::StateAndLayers BaseCkfTrajectoryBuilder::findStateAndLayers(
0228     const TrajectorySeed& seed, const TempTrajectory& traj) const {
0229   if (traj.empty()) {
0230     //set the currentState to be the one from the trajectory seed starting point
0231     PTrajectoryStateOnDet const& ptod = seed.startingState();
0232     DetId id(ptod.detId());
0233     const GeomDet* g = theMeasurementTracker->geomTracker()->idToDet(id);
0234     const Surface* surface = &g->surface();
0235
0236     TSOS currentState(
0237         trajectoryStateTransform::transientState(ptod, surface, forwardPropagator(seed)->magneticField()));
0238     const DetLayer* lastLayer = theMeasurementTracker->geometricSearchTracker()->detLayer(id);
0239     return StateAndLayers(currentState,
0240                           theNavigationSchool->nextLayers(*lastLayer, *currentState.freeState(), traj.direction()));
0241   } else {
0242     TSOS const& currentState = traj.lastMeasurement().updatedState();
0243     return StateAndLayers(
0244         currentState, theNavigationSchool->nextLayers(*traj.lastLayer(), *currentState.freeState(), traj.direction()));
0245   }
0246 }
0247
0248 void BaseCkfTrajectoryBuilder::setData(const MeasurementTrackerEvent* data) {
0249   // possibly do some sanity check here
0250   theMeasurementTracker = data;
0251 }
0252
0253 void BaseCkfTrajectoryBuilder::setEvent(const edm::Event& event) const {
0254   std::cerr << "ERROR SetEvent called on " << typeid(*this).name()
0255             << (theMeasurementTracker ? " with valid " : "witout any ") << "MeasurementTrackerEvent" << std::endl;
0256 }
0257
0258 void BaseCkfTrajectoryBuilder::unset() const {
0259   std::cerr << "ERROR unSet called on " << typeid(*this).name()
0260             << (theMeasurementTracker ? " with valid " : "witout any ") << "MeasurementTrackerEvent" << std::endl;
0261 }
0262
0263 void BaseCkfTrajectoryBuilder::setEvent(const edm::Event& iEvent,
0264                                         const edm::EventSetup& iSetup,
0265                                         const MeasurementTrackerEvent* data) {
0266   theUpdator = &iSetup.getData(theUpdatorToken);
0267   thePropagatorAlong = &iSetup.getData(thePropagatorAlongToken);
0268   thePropagatorOpposite = &iSetup.getData(thePropagatorOppositeToken);
0269   theEstimator = &iSetup.getData(theEstimatorToken);
0270   theTTRHBuilder = &iSetup.getData(theRecHitBuilderToken);
0271
0272   setData(data);
0273   if (theFilter)
0274     theFilter->setEvent(iEvent, iSetup);
0275   if (theInOutFilter)
0276     theInOutFilter->setEvent(iEvent, iSetup);
0277   setEvent_(iEvent, iSetup);
0278 }