File indexing completed on 2024-04-06 12:31:38
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "TrackingTools/PatternTools/interface/TrajectorySmoother.h"
0010 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0011 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0012
0013 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0014 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0015
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/Utilities/interface/isFinite.h"
0018 #include "TrackingTools/TransientTrackingRecHit/interface/InvalidTransientRecHit.h"
0019 #include "CommonTools/Utils/interface/DynArray.h"
0020
0021 namespace {
0022
0023 struct KFFittingSmootherParam {
0024 explicit KFFittingSmootherParam(const edm::ParameterSet& conf)
0025 : theEstimateCut(conf.getParameter<double>("EstimateCut")),
0026 theMaxFractionOutliers(conf.getParameter<double>("MaxFractionOutliers")),
0027 theMaxNumberOfOutliers(conf.getParameter<int>("MaxNumberOfOutliers")),
0028 theNoOutliersBeginEnd(conf.getParameter<bool>("NoOutliersBeginEnd")),
0029 theMinDof(conf.getParameter<int>("MinDof")),
0030 theMinNumberOfHits(conf.getParameter<int>("MinNumberOfHits")),
0031 theMinNumberOfHitsHighEta(conf.getParameter<int>("MinNumberOfHitsHighEta")),
0032 theHighEtaSwitch(conf.getParameter<double>("HighEtaSwitch")),
0033 rejectTracksFlag(conf.getParameter<bool>("RejectTracks")),
0034 breakTrajWith2ConsecutiveMissing(conf.getParameter<bool>("BreakTrajWith2ConsecutiveMissing")),
0035 noInvalidHitsBeginEnd(conf.getParameter<bool>("NoInvalidHitsBeginEnd")) {}
0036
0037 double theEstimateCut;
0038
0039 float theMaxFractionOutliers;
0040 int theMaxNumberOfOutliers;
0041 bool theNoOutliersBeginEnd;
0042 int theMinDof;
0043
0044 int theMinNumberOfHits;
0045 int theMinNumberOfHitsHighEta;
0046 double theHighEtaSwitch;
0047 bool rejectTracksFlag;
0048 bool breakTrajWith2ConsecutiveMissing;
0049 bool noInvalidHitsBeginEnd;
0050 };
0051
0052 class KFFittingSmoother final : public TrajectoryFitter, private KFFittingSmootherParam {
0053 public:
0054 ~KFFittingSmoother() override {}
0055
0056 KFFittingSmoother(const TrajectoryFitter& aFitter,
0057 const TrajectorySmoother& aSmoother,
0058 const edm::ParameterSet& conf)
0059 : KFFittingSmootherParam(conf), theFitter(aFitter.clone()), theSmoother(aSmoother.clone()) {}
0060
0061 private:
0062 static void fillDescriptions(edm::ParameterSetDescription& desc) {
0063 desc.add<double>("EstimateCut", -1);
0064 desc.add<double>("MaxFractionOutliers", 0.3);
0065 desc.add<int>("MaxNumberOfOutliers", 3);
0066 desc.add<int>("MinDof", 2);
0067 desc.add<bool>("NoOutliersBeginEnd", false);
0068 desc.add<int>("MinNumberOfHits", 5);
0069 desc.add<int>("MinNumberOfHitsHighEta", 5);
0070 desc.add<double>("HighEtaSwitch", 5.0);
0071 desc.add<bool>("RejectTracks", true);
0072 desc.add<bool>("BreakTrajWith2ConsecutiveMissing", true);
0073 desc.add<bool>("NoInvalidHitsBeginEnd", true);
0074 desc.add<double>("LogPixelProbabilityCut", 0);
0075 }
0076
0077 int getNhitCutValue(const Trajectory& t,
0078 double theHighEtaSwitch,
0079 int theMinNumberOfHitsHighEta,
0080 int theMinNumberOfHits) const {
0081 double sinhTrajEta2 = std::numeric_limits<double>::max();
0082 if (!t.empty() && t.isValid()) {
0083
0084
0085
0086
0087 double pt = t.lastMeasurement().updatedState().freeTrajectoryState()->momentum().perp();
0088 double pz = t.lastMeasurement().updatedState().freeTrajectoryState()->momentum().z();
0089 sinhTrajEta2 = (pz * pz) / (pt * pt);
0090 }
0091 double myEtaSwitch = sinh(theHighEtaSwitch);
0092 const auto thisHitCut =
0093 sinhTrajEta2 > (myEtaSwitch * myEtaSwitch) ? theMinNumberOfHitsHighEta : theMinNumberOfHits;
0094 return thisHitCut;
0095 }
0096
0097 Trajectory fitOne(const Trajectory& t, fitType type) const override;
0098 Trajectory fitOne(const TrajectorySeed& aSeed,
0099 const RecHitContainer& hits,
0100 const TrajectoryStateOnSurface& firstPredTsos,
0101 fitType type) const override;
0102 Trajectory fitOne(const TrajectorySeed& aSeed, const RecHitContainer& hits, fitType type) const override;
0103
0104 const TrajectoryFitter* fitter() const { return theFitter.get(); }
0105 const TrajectorySmoother* smoother() const { return theSmoother.get(); }
0106
0107 std::unique_ptr<TrajectoryFitter> clone() const override {
0108 return std::unique_ptr<TrajectoryFitter>(new KFFittingSmoother(*theFitter, *theSmoother, *this));
0109 }
0110
0111 void setHitCloner(TkCloner const* hc) override {
0112 theFitter->setHitCloner(hc);
0113 theSmoother->setHitCloner(hc);
0114 }
0115
0116 KFFittingSmoother(const TrajectoryFitter& aFitter,
0117 const TrajectorySmoother& aSmoother,
0118 KFFittingSmootherParam const& other)
0119 : KFFittingSmootherParam(other), theFitter(aFitter.clone()), theSmoother(aSmoother.clone()) {}
0120
0121 Trajectory smoothingStep(Trajectory&& fitted) const {
0122 const auto thisHitCut = getNhitCutValue(fitted, theHighEtaSwitch, theMinNumberOfHitsHighEta, theMinNumberOfHits);
0123
0124 if (theEstimateCut > 0) {
0125
0126 while (
0127 !fitted.empty() && fitted.foundHits() >= thisHitCut &&
0128 (!fitted.lastMeasurement().recHitR().isValid() || (fitted.lastMeasurement().recHitR().det() != nullptr &&
0129 fitted.lastMeasurement().estimate() > theEstimateCut)))
0130 fitted.pop();
0131 if (fitted.foundHits() < thisHitCut)
0132 return Trajectory();
0133 }
0134 return theSmoother->trajectory(fitted);
0135 }
0136
0137 private:
0138 const std::unique_ptr<TrajectoryFitter> theFitter;
0139 const std::unique_ptr<TrajectorySmoother> theSmoother;
0140
0141
0142 static bool checkForNans(const Trajectory& theTraj);
0143
0144 friend class KFFittingSmootherESProducer;
0145 };
0146
0147
0148
0149 #ifdef VI_DEBUG
0150 #define DPRINT(x) std::cout << x << ": "
0151 #define PRINT std::cout
0152 #else
0153 #define DPRINT(x) LogTrace(x)
0154 #define PRINT LogTrace("")
0155 #endif
0156
0157 inline Trajectory KFFittingSmoother::fitOne(const Trajectory& t, fitType type) const {
0158 if (!t.isValid())
0159 return Trajectory();
0160 return smoothingStep(theFitter->fitOne(t, type));
0161 }
0162
0163 inline bool KFFittingSmoother::checkForNans(const Trajectory& theTraj) {
0164 if (edm::isNotFinite(theTraj.chiSquared()))
0165 return false;
0166 auto const& vtm = theTraj.measurements();
0167 for (auto const& tm : vtm) {
0168 if (edm::isNotFinite(tm.estimate()))
0169 return false;
0170 auto const& v = tm.updatedState().localParameters().vector();
0171 for (int i = 0; i < 5; ++i)
0172 if (edm::isNotFinite(v[i]))
0173 return false;
0174 if (!tm.updatedState().curvilinearError().posDef())
0175 return false;
0176 auto const& m = tm.updatedState().curvilinearError().matrix();
0177 for (int i = 0; i < 5; ++i)
0178 for (int j = 0; j < (i + 1); ++j)
0179 if (edm::isNotFinite(m(i, j)))
0180 return false;
0181 }
0182 return true;
0183 }
0184
0185 namespace {
0186 inline void print(const std::string& header, const TrajectoryStateOnSurface& tsos) {
0187 DPRINT("TrackFitters") << header << tsos.globalPosition().perp() << ' ' << tsos.globalPosition().z() << ' '
0188 << 1. / tsos.signedInverseMomentum() << ' ' << 1. / tsos.transverseCurvature() << ' '
0189 << tsos.globalMomentum().eta() << std::endl;
0190 }
0191 }
0192
0193 inline Trajectory KFFittingSmoother::fitOne(const TrajectorySeed& aSeed,
0194 const RecHitContainer& hits,
0195 const TrajectoryStateOnSurface& firstPredTsos,
0196 fitType type) const {
0197 LogDebug("TrackFitters") << "In KFFittingSmoother::fit";
0198
0199 print("firstPred ", firstPredTsos);
0200
0201 if (hits.empty())
0202 return Trajectory();
0203
0204 RecHitContainer myHits = hits;
0205 Trajectory tmp_first;
0206
0207
0208 Trajectory smoothed = smoothingStep(theFitter->fitOne(aSeed, myHits, firstPredTsos));
0209
0210 do {
0211 #ifdef EDM_ML_DEBUG
0212
0213 for (unsigned int j = 0; j < myHits.size(); j++) {
0214 if (myHits[j]->det())
0215 LogTrace("TrackFitters") << "hit #:" << j + 1 << " rawId=" << myHits[j]->det()->geographicalId().rawId()
0216 << " validity=" << myHits[j]->isValid();
0217 else
0218 LogTrace("TrackFitters") << "hit #:" << j + 1 << " Hit with no Det information";
0219 }
0220 #endif
0221
0222 #if defined(VI_DEBUG) || defined(EDM_ML_DEBUG)
0223 if (smoothed.isValid()) {
0224 print("first state ", smoothed.firstMeasurement().updatedState());
0225 print("last state ", smoothed.lastMeasurement().updatedState());
0226 }
0227 #endif
0228
0229 bool hasNaN = false;
0230 const auto thisHitCut =
0231 getNhitCutValue(smoothed, theHighEtaSwitch, theMinNumberOfHitsHighEta, theMinNumberOfHits);
0232
0233 if (!smoothed.isValid() || (hasNaN = !checkForNans(smoothed)) || (smoothed.foundHits() < thisHitCut)) {
0234 if (hasNaN)
0235 edm::LogWarning("TrackNaN") << "Track has NaN or the cov is not pos-definite";
0236 if (smoothed.foundHits() < thisHitCut)
0237 LogTrace("TrackFitters") << "smoothed.foundHits()<thisHitCut";
0238 DPRINT("TrackFitters") << "smoothed invalid => trajectory rejected with nhits/chi2 " << smoothed.foundHits()
0239 << '/' << smoothed.chiSquared() << "\n";
0240 if (rejectTracksFlag) {
0241 return Trajectory();
0242 } else {
0243 std::swap(smoothed, tmp_first);
0244 DPRINT("TrackFitters") << "smoothed invalid => returning orignal trajectory with nhits/chi2 "
0245 << smoothed.foundHits() << '/' << smoothed.chiSquared() << "\n";
0246 }
0247 break;
0248 }
0249 #ifdef EDM_ML_DEBUG
0250 else {
0251 LogTrace("TrackFitters") << "dump hits after smoothing";
0252 Trajectory::DataContainer meas = smoothed.measurements();
0253 for (Trajectory::DataContainer::iterator it = meas.begin(); it != meas.end(); ++it) {
0254 LogTrace("TrackFitters") << "hit #" << meas.end() - it - 1 << " validity=" << it->recHit()->isValid()
0255 << " det=" << it->recHit()->geographicalId().rawId();
0256 }
0257 }
0258 #endif
0259
0260 if (myHits.size() != smoothed.measurements().size())
0261 DPRINT("TrackFitters") << "lost hits. before/after: " << myHits.size() << '/' << smoothed.measurements().size()
0262 << "\n";
0263
0264 if (theEstimateCut <= 0)
0265 break;
0266
0267
0268
0269 auto msize = smoothed.measurements().size();
0270 declareDynArray(unsigned int, msize, bad);
0271 unsigned int nbad = 0;
0272 unsigned int ind = 0;
0273 unsigned int lastValid = smoothed.measurements().size();
0274 for (auto const& tm : smoothed.measurements()) {
0275 if (tm.estimate() > theEstimateCut &&
0276 tm.recHitR().det() != nullptr
0277 )
0278 bad[nbad++] = ind;
0279 if (ind < lastValid && tm.recHitR().det() != nullptr && tm.recHitR().isValid())
0280 lastValid = ind;
0281 ++ind;
0282 }
0283
0284 if (0 == nbad)
0285 break;
0286
0287 DPRINT("TrackFitters") << "size/found/outliers list " << smoothed.measurements().size() << '/'
0288 << smoothed.foundHits() << ' ' << nbad << ": ";
0289 for (auto i = 0U; i < nbad; ++i)
0290 PRINT << bad[i] << ',';
0291 PRINT << std::endl;
0292
0293 if (
0294
0295 int(nbad) > theMaxNumberOfOutliers || float(nbad) > theMaxFractionOutliers * float(smoothed.foundHits())) {
0296 DPRINT("TrackFitters") << "smoothed low quality => trajectory with nhits/chi2 " << smoothed.foundHits() << '/'
0297 << smoothed.chiSquared() << "\n";
0298 PRINT << "try to remove " << lastValid << std::endl;
0299 nbad = 0;
0300
0301
0302
0303
0304
0305
0306
0307
0308
0309
0310
0311 }
0312
0313
0314 bad[nbad++] = lastValid;
0315
0316
0317
0318 assert(smoothed.measurements().size() <= myHits.size());
0319
0320 myHits.resize(smoothed.measurements().size());
0321
0322 assert(smoothed.measurements().size() == myHits.size());
0323
0324 declareDynArray(Trajectory, nbad, smoothedCand);
0325
0326 auto NHits = myHits.size();
0327 float minChi2 = std::numeric_limits<float>::max();
0328
0329 auto loc = nbad;
0330 for (auto i = 0U; i < nbad; ++i) {
0331 auto j = NHits - bad[i] - 1;
0332 assert(myHits[j]->geographicalId() == smoothed.measurements()[bad[i]].recHitR().geographicalId());
0333 auto removedHit = myHits[j];
0334 myHits[j] = std::make_shared<InvalidTrackingRecHit>(*removedHit->det(), TrackingRecHit::missing);
0335 smoothedCand[i] = smoothingStep(theFitter->fitOne(aSeed, myHits, firstPredTsos));
0336 myHits[j] = removedHit;
0337 if (smoothedCand[i].isValid() && smoothedCand[i].chiSquared() < minChi2) {
0338 minChi2 = smoothedCand[i].chiSquared();
0339 loc = i;
0340 }
0341 }
0342
0343 if (loc == nbad) {
0344 DPRINT("TrackFitters") << "New trajectories all invalid"
0345 << "\n";
0346 return Trajectory();
0347 }
0348
0349 DPRINT("TrackFitters") << "outlier removed " << bad[loc] << '/' << minChi2 << " was " << smoothed.chiSquared()
0350 << "\n";
0351
0352 if (minChi2 > smoothed.chiSquared()) {
0353 DPRINT("TrackFitters") << "removing outlier makes chi2 worse " << minChi2 << '/' << smoothed.chiSquared()
0354 << "\nOri: ";
0355 for (auto const& tm : smoothed.measurements())
0356 PRINT << tm.recHitR().geographicalId() << '/' << tm.estimate() << ' ';
0357 PRINT << "\nNew: ";
0358 for (auto const& tm : smoothedCand[loc].measurements())
0359 PRINT << tm.recHitR().geographicalId() << '/' << tm.estimate() << ' ';
0360 PRINT << "\n";
0361
0362
0363
0364 }
0365
0366 std::swap(smoothed, tmp_first);
0367 myHits[NHits - bad[loc] - 1] =
0368 std::make_shared<InvalidTrackingRecHit>(*myHits[NHits - bad[loc] - 1]->det(), TrackingRecHit::missing);
0369 std::swap(smoothed, smoothedCand[loc]);
0370
0371
0372 DPRINT("TrackFitters") << "new trajectory with nhits/chi2 " << smoothed.foundHits() << '/'
0373 << smoothed.chiSquared() << "\n";
0374
0375
0376 if (breakTrajWith2ConsecutiveMissing) {
0377 unsigned int firstinvalid = myHits.size();
0378 for (unsigned int j = 0; j < myHits.size() - 1; ++j) {
0379 if (((myHits[j]->type() == TrackingRecHit::missing) && (myHits[j]->geographicalId().rawId() != 0)) &&
0380 ((myHits[j + 1]->type() == TrackingRecHit::missing) && (myHits[j + 1]->geographicalId().rawId() != 0)) &&
0381 ((myHits[j]->geographicalId().rawId() & (~3)) !=
0382 (myHits[j + 1]->geographicalId().rawId() & (~3)))
0383 ) {
0384 firstinvalid = j;
0385 DPRINT("TrackFitters") << "Found two consecutive missing hits. First invalid: " << firstinvalid << "\n";
0386 break;
0387 }
0388 }
0389
0390
0391
0392 if (firstinvalid != myHits.size()) {
0393 myHits.erase(myHits.begin() + firstinvalid, myHits.end());
0394 smoothed = smoothingStep(theFitter->fitOne(aSeed, myHits, firstPredTsos));
0395 DPRINT("TrackFitters") << "Trajectory shortened " << smoothed.foundHits() << '/' << smoothed.chiSquared()
0396 << "\n";
0397 }
0398 }
0399
0400 }
0401 while (true);
0402
0403 if (smoothed.isValid()) {
0404 if (noInvalidHitsBeginEnd && !smoothed.empty()
0405 ) {
0406
0407 if (!smoothed.empty() && !smoothed.lastMeasurement().recHitR().isValid())
0408 LogTrace("TrackFitters") << "Last measurement is invalid";
0409
0410 while (!smoothed.empty() && !smoothed.lastMeasurement().recHitR().isValid())
0411 smoothed.pop();
0412
0413
0414 if (!smoothed.empty() && !smoothed.firstMeasurement().recHitR().isValid()) {
0415 LogTrace("TrackFitters") << "First measurement is in`valid";
0416 Trajectory tmpTraj(smoothed.seed(), smoothed.direction());
0417 Trajectory::DataContainer& meas = smoothed.measurements();
0418 auto it = meas.begin();
0419 for (; it != meas.end(); ++it)
0420 if (it->recHitR().isValid())
0421 break;
0422
0423 tmpTraj.push(std::move(*it), smoothed.chiSquared());
0424
0425 for (auto itt = it + 1; itt != meas.end(); ++itt)
0426 tmpTraj.push(std::move(*itt), 0);
0427
0428 std::swap(smoothed, tmpTraj);
0429
0430 }
0431
0432 }
0433
0434 LogTrace("TrackFitters") << "end: returning smoothed trajectory with chi2=" << smoothed.chiSquared();
0435
0436
0437
0438
0439
0440
0441
0442 }
0443
0444 return smoothed;
0445 }
0446
0447 inline Trajectory KFFittingSmoother::fitOne(const TrajectorySeed& aSeed,
0448 const RecHitContainer& hits,
0449 fitType type) const {
0450 throw cms::Exception("TrackFitters",
0451 "KFFittingSmoother::fit(TrajectorySeed, <TransientTrackingRecHit>) not implemented");
0452
0453 return Trajectory();
0454 }
0455
0456 }
0457
0458 #include "FWCore/Framework/interface/ESHandle.h"
0459
0460 #include "FWCore/Framework/interface/ESProducer.h"
0461 #include "TrackingTools/TrackFitters/interface/TrajectoryFitterRecord.h"
0462
0463 namespace {
0464
0465 class KFFittingSmootherESProducer final : public edm::ESProducer {
0466 public:
0467 KFFittingSmootherESProducer(const edm::ParameterSet& p) : pset_{p} {
0468 std::string myname = p.getParameter<std::string>("ComponentName");
0469 auto cc = setWhatProduced(this, myname);
0470 fitToken_ = cc.consumes(edm::ESInputTag("", pset_.getParameter<std::string>("Fitter")));
0471 smoothToken_ = cc.consumes(edm::ESInputTag("", pset_.getParameter<std::string>("Smoother")));
0472 }
0473
0474 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0475 edm::ParameterSetDescription desc;
0476 desc.add<std::string>("ComponentName", "KFFittingSmoother");
0477 desc.add<std::string>("Fitter", "KFFitter");
0478 desc.add<std::string>("Smoother", "KFSmoother");
0479 KFFittingSmoother::fillDescriptions(desc);
0480 descriptions.add("KFFittingSmoother", desc);
0481 }
0482
0483 std::unique_ptr<TrajectoryFitter> produce(const TrajectoryFitterRecord& iRecord) {
0484 return std::make_unique<KFFittingSmoother>(iRecord.get(fitToken_), iRecord.get(smoothToken_), pset_);
0485 }
0486
0487 private:
0488 const edm::ParameterSet pset_;
0489 edm::ESGetToken<TrajectoryFitter, TrajectoryFitterRecord> fitToken_;
0490 edm::ESGetToken<TrajectorySmoother, TrajectoryFitterRecord> smoothToken_;
0491 };
0492 }
0493
0494 #include "FWCore/Framework/interface/ModuleFactory.h"
0495
0496 DEFINE_FWK_EVENTSETUP_MODULE(KFFittingSmootherESProducer);