CkfPattern/src/CkfTrajectoryBuilder.cc

0001 #include "RecoTracker/CkfPattern/interface/CkfTrajectoryBuilder.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/ParameterSet/interface/PluginDescription.h"
0005
0006 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0007
0008 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0009 #include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h"
0010 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
0011
0012 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0013 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0014 #include "TrackingTools/TrajectoryState/interface/BasicSingleTrajectoryState.h"
0015 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0016 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0017 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0018
0019 #include "RecoTracker/CkfPattern/interface/TrajCandLess.h"
0020
0021 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0022
0023 #include "RecoTracker/CkfPattern/interface/IntermediateTrajectoryCleaner.h"
0024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0025 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0026 #include "TrackingTools/TrajectoryFiltering/interface/TrajectoryFilter.h"
0027 #include "TrackingTools/TrajectoryFiltering/interface/TrajectoryFilterFactory.h"
0028
0029 using namespace std;
0030
0031 CkfTrajectoryBuilder::CkfTrajectoryBuilder(const edm::ParameterSet& conf, edm::ConsumesCollector iC)
0032     : CkfTrajectoryBuilder(conf,
0033                            iC,
0034                            BaseCkfTrajectoryBuilder::createTrajectoryFilter(
0035                                conf.getParameter<edm::ParameterSet>("trajectoryFilter"), iC)) {}
0036
0037 CkfTrajectoryBuilder::CkfTrajectoryBuilder(const edm::ParameterSet& conf,
0038                                            edm::ConsumesCollector iC,
0039                                            std::unique_ptr<TrajectoryFilter> filter)
0040     : BaseCkfTrajectoryBuilder(conf, iC, std::move(filter)) {
0041   theMaxCand = conf.getParameter<int>("maxCand");
0042   theLostHitPenalty = conf.getParameter<double>("lostHitPenalty");
0043   theFoundHitBonus = conf.getParameter<double>("foundHitBonus");
0044   theMinHitForDoubleBonus = conf.getParameter<int>("minHitForDoubleBonus");
0045   theIntermediateCleaning = conf.getParameter<bool>("intermediateCleaning");
0046   theAlwaysUseInvalidHits = conf.getParameter<bool>("alwaysUseInvalidHits");
0047 }
0048
0049 void CkfTrajectoryBuilder::fillPSetDescription(edm::ParameterSetDescription& iDesc) {
0050   BaseCkfTrajectoryBuilder::fillPSetDescription(iDesc);
0051   iDesc.add<int>("maxCand", 5);
0052   iDesc.add<double>("lostHitPenalty", 30.);
0053   iDesc.add<double>("foundHitBonus", 0.);
0054   iDesc.add<int>("minHitForDoubleBonus", 8888);
0055   iDesc.add<bool>("intermediateCleaning", true);
0056   iDesc.add<bool>("alwaysUseInvalidHits", true);
0057
0058   edm::ParameterSetDescription psdTF;
0059   psdTF.addNode(edm::PluginDescription<TrajectoryFilterFactory>("ComponentType", true));
0060   iDesc.add<edm::ParameterSetDescription>("trajectoryFilter", psdTF);
0061 }
0062
0063 void CkfTrajectoryBuilder::setEvent_(const edm::Event& event, const edm::EventSetup& iSetup) {}
0064
0065 CkfTrajectoryBuilder::TrajectoryContainer CkfTrajectoryBuilder::trajectories(const TrajectorySeed& seed) const {
0066   TrajectoryContainer result;
0067   result.reserve(5);
0068   trajectories(seed, result);
0069   return result;
0070 }
0071
0072 void CkfTrajectoryBuilder::trajectories(const TrajectorySeed& seed,
0073                                         CkfTrajectoryBuilder::TrajectoryContainer& result) const {
0074   unsigned int tmp;
0075   buildTrajectories(seed, result, tmp, nullptr);
0076 }
0077
0078 void CkfTrajectoryBuilder::buildTrajectories(const TrajectorySeed& seed,
0079                                              TrajectoryContainer& result,
0080                                              unsigned int& nCandPerSeed,
0081                                              const TrajectoryFilter*) const {
0082   if (theMeasurementTracker == nullptr) {
0083     throw cms::Exception("LogicError")
0084         << "Asking to create trajectories to an un-initialized CkfTrajectoryBuilder.\nYou have to call clone(const "
0085            "MeasurementTrackerEvent *data) and then call trajectories on it instead.\n";
0086   }
0087
0088   TempTrajectory startingTraj = createStartingTrajectory(seed);
0089
0090   nCandPerSeed = limitedCandidates(seed, startingTraj, result);
0091 }
0092
0093 unsigned int CkfTrajectoryBuilder::limitedCandidates(const TrajectorySeed& seed,
0094                                                      TempTrajectory& startingTraj,
0095                                                      TrajectoryContainer& result) const {
0096   TempTrajectoryContainer candidates;
0097   candidates.push_back(startingTraj);
0098   std::shared_ptr<const TrajectorySeed> sharedSeed(new TrajectorySeed(seed));
0099   return limitedCandidates(sharedSeed, candidates, result);
0100 }
0101
0102 unsigned int CkfTrajectoryBuilder::limitedCandidates(const std::shared_ptr<const TrajectorySeed>& sharedSeed,
0103                                                      TempTrajectoryContainer& candidates,
0104                                                      TrajectoryContainer& result) const {
0105   unsigned int nIter = 1;
0106   unsigned int nCands = 0;  // ignore startingTraj
0107   unsigned int prevNewCandSize = 0;
0108   TempTrajectoryContainer newCand;  // = TrajectoryContainer();
0109   newCand.reserve(theMaxCand);
0110
0111   auto score = [&](TempTrajectory const& a) {
0112     auto bonus = theFoundHitBonus;
0113     bonus += a.foundHits() > theMinHitForDoubleBonus ? bonus : 0;
0114     return a.chiSquared() + a.lostHits() * theLostHitPenalty - bonus * a.foundHits();
0115   };
0116
0117   auto trajCandLess = [&](TempTrajectory const& a, TempTrajectory const& b) { return score(a) < score(b); };
0118
0119   while (!candidates.empty()) {
0120     newCand.clear();
0121     bool full = false;
0122     for (auto traj = candidates.begin(); traj != candidates.end(); traj++) {
0123       std::vector<TM> meas;
0124       findCompatibleMeasurements(*sharedSeed, *traj, meas);
0125
0126       // --- method for debugging
0127       if (!analyzeMeasurementsDebugger(
0128               *traj, meas, theMeasurementTracker, forwardPropagator(*sharedSeed), theEstimator, theTTRHBuilder))
0129         return nCands;
0130       // ---
0131
0132       if (meas.empty()) {
0133         addToResult(sharedSeed, *traj, result);
0134       } else {
0135         std::vector<TM>::const_iterator last;
0136         if (theAlwaysUseInvalidHits)
0137           last = meas.end();
0138         else {
0139           if (meas.front().recHit()->isValid()) {
0140             last = find_if(meas.begin(), meas.end(), [](auto const& meas) { return !meas.recHit()->isValid(); });
0141           } else
0142             last = meas.end();
0143         }
0144
0145         for (auto itm = meas.begin(); itm != last; itm++) {
0146           TempTrajectory newTraj = *traj;
0147           updateTrajectory(newTraj, std::move(*itm));
0148
0149           if (toBeContinued(newTraj)) {
0150             if (full) {
0151               bool better = trajCandLess(newTraj, newCand.front());
0152               if (better) {
0153                 // replace worst
0154                 std::pop_heap(newCand.begin(), newCand.end(), trajCandLess);
0155                 newCand.back().swap(newTraj);
0156                 std::push_heap(newCand.begin(), newCand.end(), trajCandLess);
0157               }  // else? no need to add it just to remove it later!
0158             } else {
0159               newCand.push_back(std::move(newTraj));
0160               full = (int)newCand.size() == theMaxCand;
0161               if (full)
0162                 std::make_heap(newCand.begin(), newCand.end(), trajCandLess);
0163             }
0164           } else {
0165             addToResult(sharedSeed, newTraj, result);
0166             //// don't know yet
0167           }
0168         }
0169       }
0170
0171       // account only new candidates, i.e.
0172       // - 1 candidate -> 1 candidate, don't increase count
0173       // - 1 candidate -> 2 candidates, increase count by 1
0174       nCands += newCand.size() - prevNewCandSize;
0175       prevNewCandSize = newCand.size();
0176
0177       assert((int)newCand.size() <= theMaxCand);
0178       if (full)
0179         assert((int)newCand.size() == theMaxCand);
0180     }  // end loop on candidates
0181
0182     // no reason to sort  (no sorting in Grouped version!)
0183     if (theIntermediateCleaning)
0184       IntermediateTrajectoryCleaner::clean(newCand);
0185
0186     candidates.swap(newCand);
0187
0188     LogDebug("CkfPattern") << result.size() << " candidates after " << nIter++ << " CKF iteration: \n"
0189                            << PrintoutHelper::dumpCandidates(result) << "\n " << candidates.size()
0190                            << " running candidates are: \n"
0191                            << PrintoutHelper::dumpCandidates(candidates);
0192   }
0193   return nCands;
0194 }
0195
0196 void CkfTrajectoryBuilder::updateTrajectory(TempTrajectory& traj, TM&& tm) const {
0197   auto&& predictedState = tm.predictedState();
0198   auto&& hit = tm.recHit();
0199   if (hit->isValid()) {
0200     auto&& upState = theUpdator->update(predictedState, *hit);
0201     traj.emplace(predictedState, std::move(upState), hit, tm.estimate(), tm.layer());
0202   } else {
0203     traj.emplace(predictedState, hit, 0, tm.layer());
0204   }
0205 }
0206
0207 void CkfTrajectoryBuilder::findCompatibleMeasurements(const TrajectorySeed& seed,
0208                                                       const TempTrajectory& traj,
0209                                                       std::vector<TrajectoryMeasurement>& result) const {
0210   int invalidHits = 0;
0211   //Use findStateAndLayers which handles the hitless seed use case
0212   std::pair<TSOS, std::vector<const DetLayer*> >&& stateAndLayers = findStateAndLayers(seed, traj);
0213   if (stateAndLayers.second.empty())
0214     return;
0215
0216   auto layerBegin = stateAndLayers.second.begin();
0217   auto layerEnd = stateAndLayers.second.end();
0218   LogDebug("CkfPattern") << "looping on " << stateAndLayers.second.size() << " layers.";
0219   const Propagator* fwdPropagator = forwardPropagator(seed);
0220   for (auto il = layerBegin; il != layerEnd; il++) {
0221     LogDebug("CkfPattern") << "looping on a layer in findCompatibleMeasurements.\n last layer: " << traj.lastLayer()
0222                            << " current layer: " << (*il);
0223
0224     TSOS stateToUse = stateAndLayers.first;
0225     //Added protection before asking for the lastLayer on the trajectory
0226     if UNLIKELY (!traj.empty() && (*il) == traj.lastLayer()) {
0227       LogDebug("CkfPattern") << " self propagating in findCompatibleMeasurements.\n from: \n" << stateToUse;
0228       //self navigation case
0229       // go to a middle point first
0230       TransverseImpactPointExtrapolator middle;
0231       GlobalPoint center(0, 0, 0);
0232       stateToUse = middle.extrapolate(stateToUse, center, *fwdPropagator);
0233
0234       if (!stateToUse.isValid())
0235         continue;
0236       LogDebug("CkfPattern") << "to: " << stateToUse;
0237     }
0238
0239     LayerMeasurements layerMeasurements(theMeasurementTracker->measurementTracker(), *theMeasurementTracker);
0240     std::vector<TrajectoryMeasurement>&& tmp =
0241         layerMeasurements.measurements((**il), stateToUse, *fwdPropagator, *theEstimator);
0242
0243     if (!tmp.empty()) {
0244       if (result.empty())
0245         result.swap(tmp);
0246       else {
0247         // keep one dummy TM at the end, skip the others
0248         result.insert(
0249             result.end() - invalidHits, std::make_move_iterator(tmp.begin()), std::make_move_iterator(tmp.end()));
0250       }
0251       invalidHits++;
0252     }
0253   }
0254
0255   // sort the final result, keep dummy measurements at the end
0256   if (result.size() > 1) {
0257     std::sort(result.begin(), result.end() - invalidHits, TrajMeasLessEstim());
0258   }
0259
0260   LogDebug("CkfPattern") << "starting from:\n"
0261                          << "x: " << stateAndLayers.first.globalPosition() << "\n"
0262                          << "p: " << stateAndLayers.first.globalMomentum() << "\n"
0263                          << PrintoutHelper::dumpMeasurements(result);
0264
0265 #ifdef DEBUG_INVALID
0266   bool afterInvalid = false;
0267   for (vector<TM>::const_iterator i = result.begin(); i != result.end(); i++) {
0268     if (!i->recHit().isValid())
0269       afterInvalid = true;
0270     if (afterInvalid && i->recHit().isValid()) {
0271       edm::LogError("CkfPattern") << "CkfTrajectoryBuilder error: valid hit after invalid!";
0272     }
0273   }
0274 #endif
0275 }