File indexing completed on 2022-01-10 06:14:16
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 theIntermediateCleaning = conf.getParameter<bool>("intermediateCleaning");
0044 theAlwaysUseInvalidHits = conf.getParameter<bool>("alwaysUseInvalidHits");
0045
0046
0047
0048
0049
0050
0051
0052 }
0053
0054 void CkfTrajectoryBuilder::fillPSetDescription(edm::ParameterSetDescription& iDesc) {
0055 BaseCkfTrajectoryBuilder::fillPSetDescription(iDesc);
0056 iDesc.add<int>("maxCand", 5);
0057 iDesc.add<double>("lostHitPenalty", 30.);
0058 iDesc.add<bool>("intermediateCleaning", true);
0059 iDesc.add<bool>("alwaysUseInvalidHits", true);
0060
0061 edm::ParameterSetDescription psdTF;
0062 psdTF.addNode(edm::PluginDescription<TrajectoryFilterFactory>("ComponentType", true));
0063 iDesc.add<edm::ParameterSetDescription>("trajectoryFilter", psdTF);
0064 }
0065
0066
0067
0068
0069
0070
0071
0072
0073 void CkfTrajectoryBuilder::setEvent_(const edm::Event& event, const edm::EventSetup& iSetup) {}
0074
0075 CkfTrajectoryBuilder::TrajectoryContainer CkfTrajectoryBuilder::trajectories(const TrajectorySeed& seed) const {
0076 TrajectoryContainer result;
0077 result.reserve(5);
0078 trajectories(seed, result);
0079 return result;
0080 }
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140 void CkfTrajectoryBuilder::trajectories(const TrajectorySeed& seed,
0141 CkfTrajectoryBuilder::TrajectoryContainer& result) const {
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156 unsigned int tmp;
0157 buildTrajectories(seed, result, tmp, nullptr);
0158 }
0159
0160 TempTrajectory CkfTrajectoryBuilder::buildTrajectories(const TrajectorySeed& seed,
0161 TrajectoryContainer& result,
0162 unsigned int& nCandPerSeed,
0163 const TrajectoryFilter*) const {
0164 if (theMeasurementTracker == nullptr) {
0165 throw cms::Exception("LogicError")
0166 << "Asking to create trajectories to an un-initialized CkfTrajectoryBuilder.\nYou have to call clone(const "
0167 "MeasurementTrackerEvent *data) and then call trajectories on it instead.\n";
0168 }
0169
0170 TempTrajectory startingTraj = createStartingTrajectory(seed);
0171
0172
0173
0174 nCandPerSeed = limitedCandidates(seed, startingTraj, result);
0175
0176 return startingTraj;
0177
0178
0179
0180
0181
0182
0183
0184 }
0185
0186 unsigned int CkfTrajectoryBuilder::limitedCandidates(const TrajectorySeed& seed,
0187 TempTrajectory& startingTraj,
0188 TrajectoryContainer& result) const {
0189 TempTrajectoryContainer candidates;
0190 candidates.push_back(startingTraj);
0191 std::shared_ptr<const TrajectorySeed> sharedSeed(new TrajectorySeed(seed));
0192 return limitedCandidates(sharedSeed, candidates, result);
0193 }
0194
0195 unsigned int CkfTrajectoryBuilder::limitedCandidates(const std::shared_ptr<const TrajectorySeed>& sharedSeed,
0196 TempTrajectoryContainer& candidates,
0197 TrajectoryContainer& result) const {
0198 unsigned int nIter = 1;
0199 unsigned int nCands = 0;
0200 unsigned int prevNewCandSize = 0;
0201 TempTrajectoryContainer newCand;
0202 newCand.reserve(2 * theMaxCand);
0203
0204 auto trajCandLess = [&](TempTrajectory const& a, TempTrajectory const& b) {
0205 return (a.chiSquared() + a.lostHits() * theLostHitPenalty) < (b.chiSquared() + b.lostHits() * theLostHitPenalty);
0206 };
0207
0208 while (!candidates.empty()) {
0209 newCand.clear();
0210 for (auto traj = candidates.begin(); traj != candidates.end(); traj++) {
0211 std::vector<TM> meas;
0212 findCompatibleMeasurements(*sharedSeed, *traj, meas);
0213
0214
0215 if (!analyzeMeasurementsDebugger(
0216 *traj, meas, theMeasurementTracker, forwardPropagator(*sharedSeed), theEstimator, theTTRHBuilder))
0217 return nCands;
0218
0219
0220 if (meas.empty()) {
0221 addToResult(sharedSeed, *traj, result);
0222 } else {
0223 std::vector<TM>::const_iterator last;
0224 if (theAlwaysUseInvalidHits)
0225 last = meas.end();
0226 else {
0227 if (meas.front().recHit()->isValid()) {
0228 last = find_if(meas.begin(), meas.end(), [](auto const& meas) { return !meas.recHit()->isValid(); });
0229 } else
0230 last = meas.end();
0231 }
0232
0233 for (auto itm = meas.begin(); itm != last; itm++) {
0234 TempTrajectory newTraj = *traj;
0235 updateTrajectory(newTraj, std::move(*itm));
0236
0237 if (toBeContinued(newTraj)) {
0238 newCand.push_back(std::move(newTraj));
0239 std::push_heap(newCand.begin(), newCand.end(), trajCandLess);
0240 } else {
0241 addToResult(sharedSeed, newTraj, result);
0242
0243 }
0244 }
0245 }
0246
0247
0248
0249
0250 nCands += newCand.size() - prevNewCandSize;
0251 prevNewCandSize = newCand.size();
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277
0278
0279
0280
0281
0282
0283
0284
0285 while ((int)newCand.size() > theMaxCand) {
0286 std::pop_heap(newCand.begin(), newCand.end(), trajCandLess);
0287
0288 newCand.pop_back();
0289 }
0290
0291
0292
0293
0294
0295
0296
0297
0298
0299
0300
0301
0302 }
0303
0304 std::sort_heap(newCand.begin(), newCand.end(), trajCandLess);
0305 if (theIntermediateCleaning)
0306 IntermediateTrajectoryCleaner::clean(newCand);
0307
0308 candidates.swap(newCand);
0309
0310 LogDebug("CkfPattern") << result.size() << " candidates after " << nIter++ << " CKF iteration: \n"
0311 << PrintoutHelper::dumpCandidates(result) << "\n " << candidates.size()
0312 << " running candidates are: \n"
0313 << PrintoutHelper::dumpCandidates(candidates);
0314 }
0315 return nCands;
0316 }
0317
0318 void CkfTrajectoryBuilder::updateTrajectory(TempTrajectory& traj, TM&& tm) const {
0319 auto&& predictedState = tm.predictedState();
0320 auto&& hit = tm.recHit();
0321 if (hit->isValid()) {
0322 auto&& upState = theUpdator->update(predictedState, *hit);
0323 traj.emplace(predictedState, std::move(upState), hit, tm.estimate(), tm.layer());
0324 } else {
0325 traj.emplace(predictedState, hit, 0, tm.layer());
0326 }
0327 }
0328
0329 void CkfTrajectoryBuilder::findCompatibleMeasurements(const TrajectorySeed& seed,
0330 const TempTrajectory& traj,
0331 std::vector<TrajectoryMeasurement>& result) const {
0332 int invalidHits = 0;
0333
0334 std::pair<TSOS, std::vector<const DetLayer*> >&& stateAndLayers = findStateAndLayers(seed, traj);
0335 if (stateAndLayers.second.empty())
0336 return;
0337
0338 auto layerBegin = stateAndLayers.second.begin();
0339 auto layerEnd = stateAndLayers.second.end();
0340 LogDebug("CkfPattern") << "looping on " << stateAndLayers.second.size() << " layers.";
0341 const Propagator* fwdPropagator = forwardPropagator(seed);
0342 for (auto il = layerBegin; il != layerEnd; il++) {
0343 LogDebug("CkfPattern") << "looping on a layer in findCompatibleMeasurements.\n last layer: " << traj.lastLayer()
0344 << " current layer: " << (*il);
0345
0346 TSOS stateToUse = stateAndLayers.first;
0347
0348 if UNLIKELY (!traj.empty() && (*il) == traj.lastLayer()) {
0349 LogDebug("CkfPattern") << " self propagating in findCompatibleMeasurements.\n from: \n" << stateToUse;
0350
0351
0352 TransverseImpactPointExtrapolator middle;
0353 GlobalPoint center(0, 0, 0);
0354 stateToUse = middle.extrapolate(stateToUse, center, *fwdPropagator);
0355
0356 if (!stateToUse.isValid())
0357 continue;
0358 LogDebug("CkfPattern") << "to: " << stateToUse;
0359 }
0360
0361 LayerMeasurements layerMeasurements(theMeasurementTracker->measurementTracker(), *theMeasurementTracker);
0362 std::vector<TrajectoryMeasurement>&& tmp =
0363 layerMeasurements.measurements((**il), stateToUse, *fwdPropagator, *theEstimator);
0364
0365 if (!tmp.empty()) {
0366 if (result.empty())
0367 result.swap(tmp);
0368 else {
0369
0370 result.insert(
0371 result.end() - invalidHits, std::make_move_iterator(tmp.begin()), std::make_move_iterator(tmp.end()));
0372 }
0373 invalidHits++;
0374 }
0375 }
0376
0377
0378 if (result.size() > 1) {
0379 std::sort(result.begin(), result.end() - invalidHits, TrajMeasLessEstim());
0380 }
0381
0382 LogDebug("CkfPattern") << "starting from:\n"
0383 << "x: " << stateAndLayers.first.globalPosition() << "\n"
0384 << "p: " << stateAndLayers.first.globalMomentum() << "\n"
0385 << PrintoutHelper::dumpMeasurements(result);
0386
0387 #ifdef DEBUG_INVALID
0388 bool afterInvalid = false;
0389 for (vector<TM>::const_iterator i = result.begin(); i != result.end(); i++) {
0390 if (!i->recHit().isValid())
0391 afterInvalid = true;
0392 if (afterInvalid && i->recHit().isValid()) {
0393 edm::LogError("CkfPattern") << "CkfTrajectoryBuilder error: valid hit after invalid!";
0394 }
0395 }
0396 #endif
0397
0398
0399 }