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   if (inOut) {
0134     // if (theInOutFilter == 0) edm::LogError("CkfPattern") << "CkfTrajectoryBuilder error: trying to use dedicated filter for in-out tracking phase, when none specified";
0135     return theInOutFilter->toBeContinued(traj);
0136   } else {
0137     return theFilter->toBeContinued(traj);
0138   }
0139 }
0140
0141 bool BaseCkfTrajectoryBuilder::qualityFilter(const TempTrajectory& traj, bool inOut) const {
0142   // Called after building a trajectory is completed, to see if it is good enough
0143   // to keep.
0144   if (inOut) {
0145     // if (theInOutFilter == 0) edm::LogError("CkfPattern") << "CkfTrajectoryBuilder error: trying to use dedicated filter for in-out tracking phase, when none specified";
0146     return theInOutFilter->qualityFilter(traj);
0147   } else {
0148     return theFilter->qualityFilter(traj);
0149   }
0150 }
0151
0152 void BaseCkfTrajectoryBuilder::addToResult(std::shared_ptr<const TrajectorySeed> const& seed,
0153                                            TempTrajectory& tmptraj,
0154                                            TrajectoryContainer& result,
0155                                            bool inOut) const {
0156   // quality check
0157   if (!qualityFilter(tmptraj, inOut))
0158     return;
0159   Trajectory traj = tmptraj.toTrajectory();
0160   traj.setSharedSeed(seed);
0161   // discard latest dummy measurements
0162   while (!traj.empty() && !traj.lastMeasurement().recHit()->isValid())
0163     traj.pop();
0164   LogDebug("CkfPattern") << inOut << "=inOut option. pushing a Trajectory with: " << traj.foundHits() << " found hits. "
0165                          << traj.lostHits()
0166                          << " lost hits. Popped :" << (tmptraj.measurements().size()) - (traj.measurements().size())
0167                          << " hits.";
0168   result.push_back(std::move(traj));
0169 }
0170
0171 void BaseCkfTrajectoryBuilder::addToResult(TempTrajectory const& tmptraj,
0172                                            TempTrajectoryContainer& result,
0173                                            bool inOut) const {
0174   // quality check
0175   if (!qualityFilter(tmptraj, inOut))
0176     return;
0177   // discard latest dummy measurements
0178   TempTrajectory traj = tmptraj;
0179   while (!traj.empty() && !traj.lastMeasurement().recHit()->isValid())
0180     traj.pop();
0181   LogDebug("CkfPattern") << inOut << "=inOut option. pushing a TempTrajectory with: " << traj.foundHits()
0182                          << " found hits. " << traj.lostHits()
0183                          << " lost hits. Popped :" << (tmptraj.measurements().size()) - (traj.measurements().size())
0184                          << " hits.";
0185   result.push_back(std::move(traj));
0186 }
0187
0188 void BaseCkfTrajectoryBuilder::moveToResult(TempTrajectory&& traj, TempTrajectoryContainer& result, bool inOut) const {
0189   // quality check
0190   if (!qualityFilter(traj, inOut))
0191     return;
0192   // discard latest dummy measurements
0193   while (!traj.empty() && !traj.lastMeasurement().recHitR().isValid())
0194     traj.pop();
0195   LogDebug("CkfPattern") << inOut << "=inOut option. pushing a TempTrajectory with: " << traj.foundHits()
0196                          << " found hits. " << traj.lostHits();
0197   //            <<" lost hits. Popped :"<<(ttraj.measurements().size())-(traj.measurements().size())<<" hits.";
0198   result.push_back(std::move(traj));
0199 }
0200
0201 BaseCkfTrajectoryBuilder::StateAndLayers BaseCkfTrajectoryBuilder::findStateAndLayers(
0202     const TrajectorySeed& seed, const TempTrajectory& traj) const {
0203   if (traj.empty()) {
0204     //set the currentState to be the one from the trajectory seed starting point
0205     PTrajectoryStateOnDet const& ptod = seed.startingState();
0206     DetId id(ptod.detId());
0207     const GeomDet* g = theMeasurementTracker->geomTracker()->idToDet(id);
0208     const Surface* surface = &g->surface();
0209
0210     TSOS currentState(
0211         trajectoryStateTransform::transientState(ptod, surface, forwardPropagator(seed)->magneticField()));
0212     const DetLayer* lastLayer = theMeasurementTracker->geometricSearchTracker()->detLayer(id);
0213     return StateAndLayers(currentState,
0214                           theNavigationSchool->nextLayers(*lastLayer, *currentState.freeState(), traj.direction()));
0215   } else {
0216     TSOS const& currentState = traj.lastMeasurement().updatedState();
0217     return StateAndLayers(
0218         currentState, theNavigationSchool->nextLayers(*traj.lastLayer(), *currentState.freeState(), traj.direction()));
0219   }
0220 }
0221
0222 void BaseCkfTrajectoryBuilder::setData(const MeasurementTrackerEvent* data) {
0223   // possibly do some sanity check here
0224   theMeasurementTracker = data;
0225 }
0226
0227 void BaseCkfTrajectoryBuilder::setEvent(const edm::Event& event) const {
0228   std::cerr << "ERROR SetEvent called on " << typeid(*this).name()
0229             << (theMeasurementTracker ? " with valid " : "witout any ") << "MeasurementTrackerEvent" << std::endl;
0230 }
0231
0232 void BaseCkfTrajectoryBuilder::unset() const {
0233   std::cerr << "ERROR unSet called on " << typeid(*this).name()
0234             << (theMeasurementTracker ? " with valid " : "witout any ") << "MeasurementTrackerEvent" << std::endl;
0235 }
0236
0237 void BaseCkfTrajectoryBuilder::setEvent(const edm::Event& iEvent,
0238                                         const edm::EventSetup& iSetup,
0239                                         const MeasurementTrackerEvent* data) {
0240   theUpdator = &iSetup.getData(theUpdatorToken);
0241   thePropagatorAlong = &iSetup.getData(thePropagatorAlongToken);
0242   thePropagatorOpposite = &iSetup.getData(thePropagatorOppositeToken);
0243   theEstimator = &iSetup.getData(theEstimatorToken);
0244   theTTRHBuilder = &iSetup.getData(theRecHitBuilderToken);
0245
0246   setData(data);
0247   if (theFilter)
0248     theFilter->setEvent(iEvent, iSetup);
0249   if (theInOutFilter)
0250     theInOutFilter->setEvent(iEvent, iSetup);
0251   setEvent_(iEvent, iSetup);
0252 }