File indexing completed on 2024-12-20 03:14:09
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;
0107 unsigned int prevNewCandSize = 0;
0108 TempTrajectoryContainer newCand;
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
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
0154 std::pop_heap(newCand.begin(), newCand.end(), trajCandLess);
0155 newCand.back().swap(newTraj);
0156 std::push_heap(newCand.begin(), newCand.end(), trajCandLess);
0157 }
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
0167 }
0168 }
0169 }
0170
0171
0172
0173
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 }
0181
0182
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
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
0226 if UNLIKELY (!traj.empty() && (*il) == traj.lastLayer()) {
0227 LogDebug("CkfPattern") << " self propagating in findCompatibleMeasurements.\n from: \n" << stateToUse;
0228
0229
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
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
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 }