diff --git a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc index c2b3d8781c546..53b05224c25ca 100644 --- a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc @@ -103,9 +103,13 @@ class MkFitOutputConverter : public edm::global::EDProducer<> { const edm::ESGetToken mkFitGeomToken_; const edm::EDPutTokenT putTrackCandidateToken_; const edm::EDPutTokenT> putSeedStopInfoToken_; - const std::string ttrhBuilderName_; - const std::string propagatorAlongName_; - const std::string propagatorOppositeName_; + + const float qualityMaxInvPt_; + const float qualityMinTheta_; + const float qualityMaxRsq_; + const float qualityMaxZ_; + const float qualityMaxPosErrSq_; + const bool qualitySignPt_; }; MkFitOutputConverter::MkFitOutputConverter(edm::ParameterSet const& iConfig) @@ -124,7 +128,13 @@ MkFitOutputConverter::MkFitOutputConverter(edm::ParameterSet const& iConfig) iConfig.getParameter("ttrhBuilder"))}, mkFitGeomToken_{esConsumes()}, putTrackCandidateToken_{produces()}, - putSeedStopInfoToken_{produces>()} {} + putSeedStopInfoToken_{produces>()}, + qualityMaxInvPt_{float(iConfig.getParameter("qualityMaxInvPt"))}, + qualityMinTheta_{float(iConfig.getParameter("qualityMinTheta"))}, + qualityMaxRsq_{float(pow(iConfig.getParameter("qualityMaxR"), 2))}, + qualityMaxZ_{float(iConfig.getParameter("qualityMaxZ"))}, + qualityMaxPosErrSq_{float(pow(iConfig.getParameter("qualityMaxPosErr"), 2))}, + qualitySignPt_{iConfig.getParameter("qualitySignPt")} {} void MkFitOutputConverter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; @@ -139,6 +149,13 @@ void MkFitOutputConverter::fillDescriptions(edm::ConfigurationDescriptions& desc desc.add("propagatorAlong", edm::ESInputTag{"", "PropagatorWithMaterial"}); desc.add("propagatorOpposite", edm::ESInputTag{"", "PropagatorWithMaterialOpposite"}); + desc.add("qualityMaxInvPt", 100)->setComment("max(1/pt) for converted tracks"); + desc.add("qualityMinTheta", 0.01)->setComment("lower bound on theta (or pi-theta) for converted tracks"); + desc.add("qualityMaxR", 120)->setComment("max(R) for the state position for converted tracks"); + desc.add("qualityMaxZ", 280)->setComment("max(|Z|) for the state position for converted tracks"); + desc.add("qualityMaxPosErr", 100)->setComment("max position error for converted tracks"); + desc.add("qualitySignPt", true)->setComment("check sign of 1/pt for converted tracks"); + descriptions.addWithDefaultLabel(desc); } @@ -194,6 +211,59 @@ TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutp LogTrace("MkFitOutputConverter") << "Candidate " << candIndex << " pT " << cand.pT() << " eta " << cand.momEta() << " phi " << cand.momPhi() << " chi2 " << cand.chi2(); + // state: check for basic quality first + if (cand.state().invpT() > qualityMaxInvPt_ || (qualitySignPt_ && cand.state().invpT() < 0) || + cand.state().theta() < qualityMinTheta_ || (M_PI - cand.state().theta()) < qualityMinTheta_ || + cand.state().posRsq() > qualityMaxRsq_ || std::abs(cand.state().z()) > qualityMaxZ_ || + (cand.state().errors.At(0, 0) + cand.state().errors.At(1, 1) + cand.state().errors.At(2, 2)) > + qualityMaxPosErrSq_) { + edm::LogInfo("MkFitOutputConverter") + << "Candidate " << candIndex << " failed state quality checks" << cand.state().parameters; + continue; + } + + auto state = cand.state(); // copy because have to modify + state.convertFromCCSToGlbCurvilinear(); + const auto& param = state.parameters; + const auto& err = state.errors; + AlgebraicSymMatrix55 cov; + for (int i = 0; i < 5; ++i) { + for (int j = i; j < 5; ++j) { + cov[i][j] = err.At(i, j); + } + } + + auto fts = FreeTrajectoryState( + GlobalTrajectoryParameters( + GlobalPoint(param[0], param[1], param[2]), GlobalVector(param[3], param[4], param[5]), state.charge, &mf), + CurvilinearTrajectoryError(cov)); + if (!fts.curvilinearError().posDef()) { + edm::LogInfo("MkFitOutputConverter") + << "Curvilinear error not pos-def\n" + << fts.curvilinearError().matrix() << "\ncandidate " << candIndex << "ignored"; + continue; + } + + //Sylvester's criterion, start from the smaller submatrix size + double det = 0; + if ((!fts.curvilinearError().matrix().Sub(0, 0).Det(det)) || det < 0) { + edm::LogInfo("MkFitOutputConverter") + << "Fail pos-def check sub2.det for candidate " << candIndex << " with fts " << fts; + continue; + } else if ((!fts.curvilinearError().matrix().Sub(0, 0).Det(det)) || det < 0) { + edm::LogInfo("MkFitOutputConverter") + << "Fail pos-def check sub3.det for candidate " << candIndex << " with fts " << fts; + continue; + } else if ((!fts.curvilinearError().matrix().Sub(0, 0).Det(det)) || det < 0) { + edm::LogInfo("MkFitOutputConverter") + << "Fail pos-def check sub4.det for candidate " << candIndex << " with fts " << fts; + continue; + } else if ((!fts.curvilinearError().matrix().Det2(det)) || det < 0) { + edm::LogInfo("MkFitOutputConverter") + << "Fail pos-def check det for candidate " << candIndex << " with fts " << fts; + continue; + } + // hits edm::OwnVector recHits; // nTotalHits() gives sum of valid hits (nFoundHits()) and invalid/missing hits. @@ -271,28 +341,6 @@ TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutp const auto seedIndex = cand.label(); LogTrace("MkFitOutputConverter") << " from seed " << seedIndex << " seed hits"; - // state - auto state = cand.state(); // copy because have to modify - state.convertFromCCSToGlbCurvilinear(); - const auto& param = state.parameters; - const auto& err = state.errors; - AlgebraicSymMatrix55 cov; - for (int i = 0; i < 5; ++i) { - for (int j = i; j < 5; ++j) { - cov[i][j] = err.At(i, j); - } - } - - auto fts = FreeTrajectoryState( - GlobalTrajectoryParameters( - GlobalPoint(param[0], param[1], param[2]), GlobalVector(param[3], param[4], param[5]), state.charge, &mf), - CurvilinearTrajectoryError(cov)); - if (!fts.curvilinearError().posDef()) { - edm::LogInfo("MkFitOutputConverter") << "Curvilinear error not pos-def\n" - << fts.curvilinearError().matrix() << "\ncandidate ignored"; - continue; - } - auto tsosDet = mkFitOutput.propagatedToFirstLayer() ? convertInnermostState(fts, recHits, propagatorAlong, propagatorOpposite)