Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add quality selections in MkFitOutputConverter #37139

Merged
merged 3 commits into from
Mar 7, 2022
Merged

add quality selections in MkFitOutputConverter #37139

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
4612033
use Sylvester's rule to check for pos-definiteness of output mkfit ca…
slava77devel Feb 18, 2022
583bb4e
add quality cuts on the track state parameters and move the covarianc…
slava77devel Feb 24, 2022
4af1441
add a quality cut on theta and position error
slava77devel Feb 25, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
100 changes: 74 additions & 26 deletions RecoTracker/MkFit/plugins/MkFitOutputConverter.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,9 +103,13 @@ class MkFitOutputConverter : public edm::global::EDProducer<> {
const edm::ESGetToken<MkFitGeometry, TrackerRecoGeometryRecord> mkFitGeomToken_;
const edm::EDPutTokenT<TrackCandidateCollection> putTrackCandidateToken_;
const edm::EDPutTokenT<std::vector<SeedStopInfo>> 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)
Expand All @@ -124,7 +128,13 @@ MkFitOutputConverter::MkFitOutputConverter(edm::ParameterSet const& iConfig)
iConfig.getParameter<edm::ESInputTag>("ttrhBuilder"))},
mkFitGeomToken_{esConsumes<MkFitGeometry, TrackerRecoGeometryRecord>()},
putTrackCandidateToken_{produces<TrackCandidateCollection>()},
putSeedStopInfoToken_{produces<std::vector<SeedStopInfo>>()} {}
putSeedStopInfoToken_{produces<std::vector<SeedStopInfo>>()},
qualityMaxInvPt_{float(iConfig.getParameter<double>("qualityMaxInvPt"))},
qualityMinTheta_{float(iConfig.getParameter<double>("qualityMinTheta"))},
qualityMaxRsq_{float(pow(iConfig.getParameter<double>("qualityMaxR"), 2))},
qualityMaxZ_{float(iConfig.getParameter<double>("qualityMaxZ"))},
qualityMaxPosErrSq_{float(pow(iConfig.getParameter<double>("qualityMaxPosErr"), 2))},
qualitySignPt_{iConfig.getParameter<bool>("qualitySignPt")} {}

void MkFitOutputConverter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
Expand All @@ -139,6 +149,13 @@ void MkFitOutputConverter::fillDescriptions(edm::ConfigurationDescriptions& desc
desc.add("propagatorAlong", edm::ESInputTag{"", "PropagatorWithMaterial"});
desc.add("propagatorOpposite", edm::ESInputTag{"", "PropagatorWithMaterialOpposite"});

desc.add<double>("qualityMaxInvPt", 100)->setComment("max(1/pt) for converted tracks");
desc.add<double>("qualityMinTheta", 0.01)->setComment("lower bound on theta (or pi-theta) for converted tracks");
desc.add<double>("qualityMaxR", 120)->setComment("max(R) for the state position for converted tracks");
desc.add<double>("qualityMaxZ", 280)->setComment("max(|Z|) for the state position for converted tracks");
desc.add<double>("qualityMaxPosErr", 100)->setComment("max position error for converted tracks");
desc.add<bool>("qualitySignPt", true)->setComment("check sign of 1/pt for converted tracks");

descriptions.addWithDefaultLabel(desc);
}

Expand Down Expand Up @@ -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<AlgebraicSymMatrix22>(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<AlgebraicSymMatrix33>(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<AlgebraicSymMatrix44>(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<TrackingRecHit> recHits;
// nTotalHits() gives sum of valid hits (nFoundHits()) and invalid/missing hits.
Expand Down Expand Up @@ -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)
Expand Down