Skip to content

Commit

Permalink
feat(gx2f): process and count holes + unit test (#3112)
Browse files Browse the repository at this point in the history
Adds the consideration of holes to the GX2F. Later, we also want to have the material in the same place, where we consider the holes.

## physmon changes
- nHoles_vs_eta
- nHoles_vs_pT
- nStates_vs_eta
- nStates_vs_pT

They make sense to me, since we are adding holes and more states.

## Blocked by:
- #3100
- #3118
- #3119
  • Loading branch information
AJPfleger authored Apr 23, 2024
1 parent d1511d0 commit 45e76e7
Show file tree
Hide file tree
Showing 3 changed files with 221 additions and 6 deletions.
Binary file modified CI/physmon/reference/performance_gx2f.root
Binary file not shown.
116 changes: 110 additions & 6 deletions Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -399,6 +399,11 @@ class Gx2Fitter {

// End the propagation and return to the fitter
if (result.finished) {
// Remove the missing surfaces that occur after the last measurement
if (result.measurementStates > 0) {
result.missedActiveSurfaces.resize(result.measurementHoles);
}

return;
}

Expand All @@ -415,15 +420,13 @@ class Gx2Fitter {
// Update:
// - Waiting for a current surface
auto surface = navigator.currentSurface(state.navigation);
if (surface != nullptr &&
surface->associatedDetectorElement() != nullptr) {
if (surface != nullptr) {
++result.surfaceCount;
ACTS_VERBOSE("Surface " << surface->geometryId() << " detected.");

// check if measurement surface
auto sourcelink_it = inputMeasurements->find(surface->geometryId());

if (sourcelink_it != inputMeasurements->end()) {
// Check if we have a measurement surface
if (auto sourcelink_it = inputMeasurements->find(surface->geometryId());
sourcelink_it != inputMeasurements->end()) {
ACTS_VERBOSE("Measurement surface " << surface->geometryId()
<< " detected.");

Expand Down Expand Up @@ -528,6 +531,107 @@ class Gx2Fitter {

// We count the processed state
++result.processedStates;

// Update the number of holes count only when encountering a
// measurement
result.measurementHoles = result.missedActiveSurfaces.size();
} else if (surface->associatedDetectorElement() != nullptr ||
surface->surfaceMaterial() != nullptr) {
// Here we handle material and holes
// TODO add material handling
ACTS_VERBOSE("Non-Measurement surface " << surface->geometryId()
<< " detected.");

// We only create track states here if there is already a measurement
// detected or if the surface has material (no holes before the first
// measurement)
if (result.measurementStates > 0
// || surface->surfaceMaterial() != nullptr
) {
ACTS_VERBOSE("Handle hole.");

auto& fittedStates = *result.fittedStates;

// Mask for the track states. We don't need Smoothed and Filtered
TrackStatePropMask mask = TrackStatePropMask::Predicted |
TrackStatePropMask::Jacobian |
TrackStatePropMask::Calibrated;

ACTS_VERBOSE(" processSurface: addTrackState");

// Add a <mask> TrackState entry multi trajectory. This allocates
// storage for all components, which we will set later.
typename traj_t::TrackStateProxy trackStateProxy =
fittedStates.makeTrackState(mask, result.lastTrackIndex);
std::size_t currentTrackIndex = trackStateProxy.index();
{
// Set the trackStateProxy components with the state from the
// ongoing propagation
{
trackStateProxy.setReferenceSurface(surface->getSharedPtr());
// Bind the transported state to the current surface
auto res = stepper.boundState(state.stepping, *surface, false,
freeToBoundCorrection);
if (!res.ok()) {
result.result = res.error();
return;
}
const auto& [boundParams, jacobian, pathLength] = *res;

// Fill the track state
trackStateProxy.predicted() = boundParams.parameters();
trackStateProxy.predictedCovariance() = state.stepping.cov;

trackStateProxy.jacobian() = jacobian;
trackStateProxy.pathLength() = pathLength;
}

// Get and set the type flags
auto typeFlags = trackStateProxy.typeFlags();
typeFlags.set(TrackStateFlag::ParameterFlag);
if (surface->surfaceMaterial() != nullptr) {
typeFlags.set(TrackStateFlag::MaterialFlag);
}

// Set hole only, if we are on a sensitive surface
if (surface->associatedDetectorElement() != nullptr) {
ACTS_VERBOSE("Detected hole on " << surface->geometryId());
// If the surface is sensitive, set the hole type flag
typeFlags.set(TrackStateFlag::HoleFlag);
} else {
ACTS_VERBOSE("Detected in-sensitive surface "
<< surface->geometryId());
}
}

ACTS_VERBOSE(
"Actor - indices after processing, before over writing:"
<< "\n "
<< "result.lastMeasurementIndex: "
<< result.lastMeasurementIndex << "\n "
<< "trackStateProxy.index(): " << trackStateProxy.index()
<< "\n "
<< "result.lastTrackIndex: " << result.lastTrackIndex
<< "\n "
<< "currentTrackIndex: " << currentTrackIndex)
result.lastTrackIndex = currentTrackIndex;

if (trackStateProxy.typeFlags().test(TrackStateFlag::HoleFlag)) {
// Count the missed surface
result.missedActiveSurfaces.push_back(surface);
}

++result.processedStates;
} else {
ACTS_VERBOSE("Ignoring hole, because no preceding measurements.");
}

if (surface->surfaceMaterial() != nullptr) {
// TODO write similar to KF?
// Update state and stepper with material effects
// materialInteractor(surface, state, stepper, navigator,
// MaterialUpdateStage::FullUpdate);
}
} else {
ACTS_INFO("Actor: This case is not implemented yet")
}
Expand Down
111 changes: 111 additions & 0 deletions Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -802,5 +802,116 @@ BOOST_AUTO_TEST_CASE(NotEnoughMeasurements) {

ACTS_INFO("*** Test: NotEnoughMeasurements -- Finish");
}

BOOST_AUTO_TEST_CASE(FindHoles) {
ACTS_INFO("*** Test: FindHoles -- Start");

std::default_random_engine rng(42);

ACTS_DEBUG("Create the detector");
// const std::size_t nSurfaces = 7;
const std::size_t nSurfaces = 8;
Detector detector;
detector.geometry = makeToyDetector(geoCtx, nSurfaces);

ACTS_DEBUG("Set the start parameters for measurement creation and fit");
const auto parametersMeasurements = makeParameters();
const auto startParametersFit = makeParameters(
7_mm, 11_mm, 15_mm, 42_ns, 10_degree, 80_degree, 1_GeV, 1_e);

ACTS_DEBUG("Create the measurements");
using SimPropagator =
Acts::Propagator<Acts::StraightLineStepper, Acts::Navigator>;
const SimPropagator simPropagator = makeStraightPropagator(detector.geometry);
const auto measurements =
createMeasurements(simPropagator, geoCtx, magCtx, parametersMeasurements,
resMapAllPixel, rng);
auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
ACTS_VERBOSE("sourceLinks.size() [before] = " << sourceLinks.size());

// We remove the first measurement in the list. This does not create a hole.
sourceLinks.erase(std::next(sourceLinks.begin(), 0));
ACTS_VERBOSE(
"sourceLinks.size() [after first erase] = " << sourceLinks.size());

// We remove the last measurement in the list. This does not create a hole.
sourceLinks.pop_back();
ACTS_VERBOSE("sourceLinks.size() [after pop] = " << sourceLinks.size());

// We remove the second to last measurement in the list. This effectivly
// creates a hole on that surface.
const std::size_t indexHole = sourceLinks.size() - 2;
ACTS_VERBOSE("Remove measurement " << indexHole);
sourceLinks.erase(std::next(sourceLinks.begin(), indexHole));
ACTS_VERBOSE("sourceLinks.size() [after second-to-last erase]= "
<< sourceLinks.size());

// We removed 3 measurements
// const std::size_t nMeasurements = nSurfaces - 2;
const std::size_t nMeasurements = nSurfaces - 3;
BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);

ACTS_DEBUG("Set up the fitter");
const Surface* rSurface = &parametersMeasurements.referenceSurface();

using RecoStepper = EigenStepper<>;
const auto recoPropagator =
makeConstantFieldPropagator<RecoStepper>(detector.geometry, 0_T);

using RecoPropagator = decltype(recoPropagator);
using Gx2Fitter =
Experimental::Gx2Fitter<RecoPropagator, VectorMultiTrajectory>;
const Gx2Fitter fitter(recoPropagator, gx2fLogger->clone());

Experimental::Gx2FitterExtensions<VectorMultiTrajectory> extensions;
extensions.calibrator
.connect<&testSourceLinkCalibrator<VectorMultiTrajectory>>();
TestSourceLink::SurfaceAccessor surfaceAccessor{*detector.geometry};
extensions.surfaceAccessor
.connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor);

const Experimental::Gx2FitterOptions gx2fOptions(
geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface,
false, false, FreeToBoundCorrection(false), 20, true, 1e-5);

Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
Acts::VectorMultiTrajectory{}};

ACTS_DEBUG("Fit the track");
ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements);
ACTS_VERBOSE("startParameter fit:\n" << startParametersFit);
const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(),
startParametersFit, gx2fOptions, tracks);

BOOST_REQUIRE(res.ok());

const auto& track = *res;

// -1, because the index starts at 0
// -2, because the first and the last surface are not part of the track
BOOST_CHECK_EQUAL(track.tipIndex(), nSurfaces - 1 - 2);
BOOST_CHECK(track.hasReferenceSurface());

// Track quantities
CHECK_CLOSE_ABS(track.chi2(), 6.5, 2.);
BOOST_CHECK_EQUAL(track.nDoF(), 10u);
BOOST_CHECK_EQUAL(track.nHoles(), 1u);
BOOST_CHECK_EQUAL(track.nMeasurements(), nMeasurements);
BOOST_CHECK_EQUAL(track.nSharedHits(), 0u);
BOOST_CHECK_EQUAL(track.nOutliers(), 0u);

// Parameters
// We need quite coarse checks here, since on different builds
// the created measurements differ in the randomness
BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -11., 7e0);
BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc1], -15., 6e0);
BOOST_CHECK_CLOSE(track.parameters()[eBoundPhi], 1e-5, 1e3);
BOOST_CHECK_CLOSE(track.parameters()[eBoundTheta], M_PI / 2, 1e-3);
BOOST_CHECK_EQUAL(track.parameters()[eBoundQOverP], 1);
BOOST_CHECK_CLOSE(track.parameters()[eBoundTime], 12591.2832360000, 1e-6);
BOOST_CHECK_CLOSE(track.covariance().determinant(), 4.7e-28, 2e0);

ACTS_INFO("*** Test: FindHoles -- Finish");
}
BOOST_AUTO_TEST_SUITE_END()
} // namespace Acts::Test

0 comments on commit 45e76e7

Please sign in to comment.