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

feat: Marking Edgeholes during track finding #3988

Draft
wants to merge 6 commits into
base: main
Choose a base branch
from
+121 −7
Draft

feat: Marking Edgeholes during track finding #3988

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
4f0fe37
Add boundary tolerance for intersection tests and inverse tolerance f…
Nov 25, 2024
ad0e897
Restore no check on boundaries in navigator
pbutti Dec 3, 2024
e9a20e6
Merge branch 'acts-project:main' into fix-edgeholes
pbutti Dec 5, 2024
5cb0ada
Merge branch 'acts-project:main' into fix-edgeholes
pbutti Dec 5, 2024
29ed187
Merge branch 'acts-project:main' into fix-edgeholes
pbutti Dec 9, 2024
fbc3b2f
Support for skipping edge holes
pbutti Dec 13, 2024
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
11 changes: 9 additions & 2 deletions Core/include/Acts/EventData/MultiTrajectoryHelpers.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ struct TrajectoryState {
std::size_t nMeasurements = 0;
std::size_t nOutliers = 0;
std::size_t nHoles = 0;
//std::size_t nEdgeHoles = 0;
double chi2Sum = 0;
std::vector<double> measurementChi2 = {};
std::vector<double> outlierChi2 = {};
Expand Down Expand Up @@ -61,7 +62,10 @@ TrajectoryState trajectoryState(const traj_t& multiTraj, std::size_t tipIndex) {
auto typeFlags = state.typeFlags();
if (typeFlags.test(Acts::TrackStateFlag::HoleFlag)) {
trajState.nHoles++;
} else if (typeFlags.test(Acts::TrackStateFlag::OutlierFlag)) {
} //else if (typeFlags.test(Acts::TrackStateFlag::EdgeHoleFlag)) {
//trajState.nEdgeHoles++;
//}
else if (typeFlags.test(Acts::TrackStateFlag::OutlierFlag)) {
trajState.nOutliers++;
trajState.outlierChi2.push_back(state.chi2());
trajState.outlierVolume.push_back(volume);
Expand Down Expand Up @@ -113,7 +117,10 @@ VolumeTrajectoryStateContainer trajectoryState(
auto typeFlags = state.typeFlags();
if (typeFlags.test(Acts::TrackStateFlag::HoleFlag)) {
trajState.nHoles++;
} else if (typeFlags.test(Acts::TrackStateFlag::OutlierFlag)) {
} //else if (typeFlags.test(Acts::TrackStateFlag::EdgeHoleFlag)) {
//trajState.nEdgeHoles++;
//}
else if (typeFlags.test(Acts::TrackStateFlag::OutlierFlag)) {
trajState.nOutliers++;
trajState.outlierChi2.push_back(state.chi2());
trajState.outlierVolume.push_back(volume);
Expand Down
21 changes: 19 additions & 2 deletions Core/include/Acts/EventData/TrackProxy.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -349,12 +349,28 @@ class TrackProxy {
return component<unsigned int, hashString("nHoles")>();
}

/// Return the number of measurements for the track. Const version
/// @return The number of measurements
/// Return the number of holes for the track. Const version
/// @return The number of holes
unsigned int nHoles() const {
return component<unsigned int, hashString("nHoles")>();
}

/// Return a mutable reference to the number of edge holes for the track.
/// Mutable version
/// @note Only available if the track proxy is not read-only
/// @return The number of edge holes
//unsigned int& nEdgeHoles()
// requires(!ReadOnly)
//{
// return component<unsigned int, hashString("nEdgeHoles")>();
//}

/// Return the number of edge holes for the track. Const version
/// @return The number of edge holes
//unsigned int nEdgeHoles() const {
// return component<unsigned int, hashString("nEdgeHoles")>();
//}

/// Return a mutable reference to the number of outliers for the track.
/// Mutable version
/// @note Only available if the track proxy is not read-only
Expand Down Expand Up @@ -596,6 +612,7 @@ class TrackProxy {

nMeasurements() = other.nMeasurements();
nHoles() = other.nHoles();
//nEdgeHoles() = other.nEdgeHoles();
nOutliers() = other.nOutliers();
nSharedHits() = other.nSharedHits();
chi2() = other.chi2();
Expand Down
3 changes: 2 additions & 1 deletion Core/include/Acts/EventData/TrackStateType.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,8 @@ enum TrackStateFlag {
HoleFlag = 3,
MaterialFlag = 4,
SharedHitFlag = 5,
NumTrackStateFlags = 6
EdgeHoleFlag = 6,
NumTrackStateFlags = 7
};

class ConstTrackStateType;
Expand Down
5 changes: 5 additions & 0 deletions Core/include/Acts/EventData/VectorTrackContainer.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,8 @@ class VectorTrackContainerBase {
return &instance.m_nMeasurements[itrack];
case "nHoles"_hash:
return &instance.m_nHoles[itrack];
//case "nEdgeHoles"_hash:
// return &instance.m_nEdgeHoles[itrack];
case "chi2"_hash:
return &instance.m_chi2[itrack];
case "ndf"_hash:
Expand Down Expand Up @@ -126,6 +128,8 @@ class VectorTrackContainerBase {
assert(result);
result = result && m_nHoles.size() == size;
assert(result);
//result = result && m_nEdgeHoles.size() == size;
//assert(result);
result = result && m_chi2.size() == size;
assert(result);
result = result && m_ndf.size() == size;
Expand Down Expand Up @@ -180,6 +184,7 @@ class VectorTrackContainerBase {

std::vector<unsigned int> m_nMeasurements;
std::vector<unsigned int> m_nHoles;
//std::vector<unsigned int> m_nEdgeHoles;
std::vector<float> m_chi2;
std::vector<unsigned int> m_ndf;
std::vector<unsigned int> m_nOutliers;
Expand Down
3 changes: 3 additions & 0 deletions Core/include/Acts/Propagator/PropagatorOptions.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,9 @@ struct PurePropagatorPlainOptions {
/// @note ignored if empty
/// @note requires `VolumeConstraintAborter` aborter
std::vector<std::uint32_t> endOfWorldVolumeIds;
/// Placeholder:: edgeHoles
bool keepEdgeHoles=false;

};

} // namespace detail
Expand Down
7 changes: 6 additions & 1 deletion Core/include/Acts/Surfaces/AnnulusBounds.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -195,7 +195,7 @@ class AnnulusBounds : public DiscBounds {
/// @return boolean indicator for the success of this operation
virtual bool inside(const Vector2& lposition, double tolR,
double tolPhi) const final;

/// Transform the strip cartesian
/// into the module polar system
///
Expand All @@ -209,6 +209,11 @@ class AnnulusBounds : public DiscBounds {

/// Private helper method
double squaredNorm(const Vector2& v, const SquareMatrix2& weight) const;

private:
bool inside(const Vector2& lposition, double tolR,
double tolPhiR,
double scale) const;
Comment on lines +214 to +216
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this can be replaced by a proper implementation of the inside check given a tolerance

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

passing the "scale" parameter avoids a division. But anyway the way it is implemented duplicates code. So, as mentioned at a different place. This was the easiest way to implement it and was anyway more meant to be a proof of concept. And so far we only see that the concept does not work all that great ;-)

};

inline SurfaceBounds::BoundsType AnnulusBounds::type() const {
Expand Down
25 changes: 24 additions & 1 deletion Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -518,6 +518,10 @@ class CombinatorialKalmanFilter {
/// Calibration context for the finding run
const CalibrationContext* calibrationContextPtr{nullptr};


/// PlaceHolder for EdgeHoles
bool keepEdgeHoles = false;

/// @brief CombinatorialKalmanFilter actor operation
///
/// @tparam propagator_state_t Type of the Propagator state
Expand Down Expand Up @@ -810,6 +814,7 @@ class CombinatorialKalmanFilter {
currentBranch = result.activeBranches.back();
prevTip = currentBranch.tipIndex();
} else {

if (expectMeasurements) {
ACTS_VERBOSE("Detected hole after measurement selection on surface "
<< surface->geometryId());
Expand All @@ -823,7 +828,23 @@ class CombinatorialKalmanFilter {
currentBranch.tipIndex() = currentTip;
auto currentState = currentBranch.outermostTrackState();
if (expectMeasurements) {
currentBranch.nHoles()++;
static const BoundaryTolerance exclude_sensor_border
= BoundaryTolerance(BoundaryTolerance::AbsoluteCartesian{-2*UnitConstants::mm,-2*UnitConstants::mm});
if (currentState.referenceSurface().insideBounds(currentState.parameters().template head<2>(),
exclude_sensor_border)) {
currentBranch.nHoles()++;
}
else {
if (!keepEdgeHoles) {
currentBranch.nHoles()++;
}
else {
auto typeFlags = currentState.typeFlags();
typeFlags.reset(TrackStateFlag::HoleFlag);
typeFlags.set(TrackStateFlag::EdgeHoleFlag);
//currentBranch.nEdgeHoles()++;
}
}
Comment on lines +831 to +847
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

IMO this is an experiment choice and it should not be part of the Core CKF. This can be put into the track state creator to put the user in control

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

IMO this is an experiment choice and it should not be part of the Core CKF. This can be put into the track state creator to put the user in control

Not sure whether there is a universally accepted definition of what a hole is. Is it really a how if the edge is within the uncertainties of the trajectory ? It is a definition, but I guess one can argue either way, whether it should be counted as a hole or not. Sure an ad hoc 2mm tolerance is not universal.

}

BranchStopperResult branchStopperResult =
Expand Down Expand Up @@ -1198,6 +1219,8 @@ class CombinatorialKalmanFilter {
combKalmanActor.updaterLogger = m_updaterLogger.get();
combKalmanActor.calibrationContextPtr = &tfOptions.calibrationContext.get();

combKalmanActor.keepEdgeHoles = tfOptions.propagatorPlainOptions.keepEdgeHoles;

// copy source link accessor, calibrator and measurement selector
combKalmanActor.m_sourceLinkAccessor = tfOptions.sourceLinkAccessor;
combKalmanActor.m_extensions = tfOptions.extensions;
Expand Down
44 changes: 44 additions & 0 deletions Core/src/Surfaces/AnnulusBounds.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -171,6 +171,44 @@ bool Acts::AnnulusBounds::inside(const Vector2& lposition, double tolR,
return true;
}

bool Acts::AnnulusBounds::inside(const Vector2& lposition, double tolR,
double tolPhiR,
double scale) const {
// locpo is PC in STRIP SYSTEM
// need to perform internal rotation induced by average phi
Vector2 locpo_rotated = m_rotationStripPC * lposition;
double phiLoc = locpo_rotated[eBoundLoc1];
double phiLocR = phiLoc*scale;
double rLoc = locpo_rotated[eBoundLoc0];

if (phiLocR < (get(eMinPhiRel)*scale - tolPhiR) ||
phiLocR > (get(eMaxPhiRel)*scale + tolPhiR)) {
return false;
}

// calculate R in MODULE SYSTEM to evaluate R-bounds
if (tolR == 0.) {
// don't need R, can use R^2
double r_mod2 =
m_shiftPC[eBoundLoc0] * m_shiftPC[eBoundLoc0] + rLoc * rLoc +
2 * m_shiftPC[eBoundLoc0] * rLoc * cos(phiLoc - m_shiftPC[eBoundLoc1]);

if (r_mod2 < get(eMinR) * get(eMinR) || r_mod2 > get(eMaxR) * get(eMaxR)) {
return false;
}
} else {
// use R
double r_mod = sqrt(
m_shiftPC[eBoundLoc0] * m_shiftPC[eBoundLoc0] + rLoc * rLoc +
2 * m_shiftPC[eBoundLoc0] * rLoc * cos(phiLoc - m_shiftPC[eBoundLoc1]));

if (r_mod < (get(eMinR) - tolR) || r_mod > (get(eMaxR) + tolR)) {
return false;
}
}
return true;
}

bool Acts::AnnulusBounds::inside(
const Vector2& lposition,
const BoundaryTolerance& boundaryTolerance) const {
Expand All @@ -181,6 +219,12 @@ bool Acts::AnnulusBounds::inside(
if (boundaryTolerance.isNone()) {
return inside(lposition, 0., 0.);
}
if (boundaryTolerance.hasAbsoluteCartesian()) {
const double r = lposition(0,0);
const BoundaryTolerance::AbsoluteCartesian& tolerance = boundaryTolerance.asAbsoluteCartesian();
return inside(lposition, tolerance.tolerance0,
tolerance.tolerance1,r);
}

if (auto absoluteBound = boundaryTolerance.asAbsoluteBoundOpt();
absoluteBound.has_value()) {
Expand Down
9 changes: 9 additions & 0 deletions Core/src/Surfaces/RectangleBounds.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,15 @@
bool Acts::RectangleBounds::inside(
const Acts::Vector2& lposition,
const Acts::BoundaryTolerance& boundaryTolerance) const {
if (boundaryTolerance.hasAbsoluteCartesian()) {

const BoundaryTolerance::AbsoluteCartesian& tolerance = boundaryTolerance.asAbsoluteCartesian();
return lposition(0,0) > (m_min.x()-tolerance.tolerance0)
&& lposition(0,0) < (m_max.x()+tolerance.tolerance0)
&& lposition(1,0) > (m_min.y()-tolerance.tolerance1)
&& lposition(1,0) < (m_max.y()+tolerance.tolerance1);
}

Comment on lines +19 to +27
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why was this necessary?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@goetzgaycken Can you comment on this?

Copy link
Contributor

@goetzgaycken goetzgaycken Dec 16, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It was just easier to implement the negative tolerance, this way then trying to add it to the more generic implementation.

Btw. the implementation was not meant for a PR, it was more a proof of concept.

Copy link
Member

@paulgessinger paulgessinger Dec 16, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doesn't insideAlignedBox do exactly the same thing in a more general way?

Copy link
Contributor

@goetzgaycken goetzgaycken Dec 18, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  • I have to admit that I do not understand how the method handles the case that the point is inside the polygon. Will computeClosestPointOnPolygon return the point itself if it is inside the polygon ? If not than I would not understand how the tolerance works .
  • anyway currently tolerance.isTolerated() does not handle a negative tolerance, and the passed distance does not carry any information whether the point is inside or outside unless "inside" will always yield a zero distance, but then it would not be possible to handle a negative tolerance.
  • far more complicated than this simple set of 4 comparisons.

return detail::insideAlignedBox(m_min, m_max, boundaryTolerance, lposition,
std::nullopt);
}
Expand Down
Loading