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

Jump to bottom

refactor(gx2f): analyse temporary track in an extra function #3799

Merged
merged 2 commits into from
Oct 31, 2024
Merged

refactor(gx2f): analyse temporary track in an extra function #3799

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
49693bd
make function
AJPfleger Oct 31, 2024
8b9b73c
Merge branch 'main' into gx2f-ts-loop
AJPfleger Oct 31, 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
174 changes: 101 additions & 73 deletions Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include "Acts/EventData/SourceLink.hpp"
#include "Acts/EventData/TrackContainerFrontendConcept.hpp"
#include "Acts/EventData/TrackParameters.hpp"
#include "Acts/EventData/TrackProxyConcept.hpp"
#include "Acts/EventData/VectorMultiTrajectory.hpp"
#include "Acts/EventData/VectorTrackContainer.hpp"
#include "Acts/Geometry/GeometryContext.hpp"
Expand Down Expand Up @@ -485,6 +486,103 @@ void addMaterialToGx2fSums(
return;
}

/// @brief Fill the GX2F system with data from a track
///
/// This function processes a track proxy and updates the aMatrix, bVector, and
/// chi2 values for the GX2F fitting system. It considers material only if
/// multiple scattering is enabled.
///
/// @tparam track_proxy_t The type of the track proxy
///
/// @param track A mutable track proxy to operate on
/// @param aMatrixExtended The aMatrix, summing over second derivatives for the fitting system
/// @param bVectorExtended The bVector, summing over first derivatives for the fitting system
/// @param chi2sum Accumulated chi2 value of the system
/// @param countNdf The number of degrees of freedom counted so far
/// @param multipleScattering Flag to consider multiple scattering in the calculation
/// @param scatteringMap Map of geometry identifiers to scattering properties, containing all scattering angles and covariances
/// @param geoIdVector A vector to store geometry identifiers for tracking processed elements
/// @param logger A logger instance
template <TrackProxyConcept track_proxy_t>
void fillGx2fSystem(
const track_proxy_t track, Eigen::MatrixXd& aMatrixExtended,
Eigen::VectorXd& bVectorExtended, double& chi2sum, std::size_t& countNdf,
const bool multipleScattering,
const std::unordered_map<GeometryIdentifier, ScatteringProperties>&
scatteringMap,
std::vector<GeometryIdentifier>& geoIdVector, const Logger& logger) {
std::vector<BoundMatrix> jacobianFromStart;
jacobianFromStart.emplace_back(BoundMatrix::Identity());

for (const auto& trackState : track.trackStates()) {
// Get and store geoId for the current surface
const GeometryIdentifier geoId = trackState.referenceSurface().geometryId();
ACTS_DEBUG("Start to investigate trackState on surface " << geoId);
const auto typeFlags = trackState.typeFlags();
const bool stateHasMeasurement =
typeFlags.test(TrackStateFlag::MeasurementFlag);
const bool stateHasMaterial = typeFlags.test(TrackStateFlag::MaterialFlag);

// First we figure out, if we would need to look into material
// surfaces at all. Later, we also check, if the material slab is
// valid, otherwise we modify this flag to ignore the material
// completely.
bool doMaterial = multipleScattering && stateHasMaterial;
if (doMaterial) {
const auto scatteringMapId = scatteringMap.find(geoId);
assert(scatteringMapId != scatteringMap.end() &&
"No scattering angles found for material surface.");
doMaterial = doMaterial && scatteringMapId->second.materialIsValid();
}

// We only consider states with a measurement (and/or material)
if (!stateHasMeasurement && !doMaterial) {
ACTS_DEBUG(" Skip state.");
continue;
}

// update all Jacobians from start
for (auto& jac : jacobianFromStart) {
jac = trackState.jacobian() * jac;
}

// Handle measurement
if (stateHasMeasurement) {
ACTS_DEBUG(" Handle measurement.");

const auto measDim = trackState.calibratedSize();

if (measDim < 1 || 6 < measDim) {
ACTS_ERROR("Can not process state with measurement with "
<< measDim << " dimensions.");
throw std::domain_error(
"Found measurement with less than 1 or more than 6 dimension(s).");
}

countNdf += measDim;

visit_measurement(measDim, [&](auto N) {
addMeasurementToGx2fSums<N>(aMatrixExtended, bVectorExtended, chi2sum,
jacobianFromStart, trackState, logger);
});
}

// Handle material
if (doMaterial) {
ACTS_DEBUG(" Handle material");
// Add for this material a new Jacobian, starting from this surface.
jacobianFromStart.emplace_back(BoundMatrix::Identity());

// Add the material contribution to the system
addMaterialToGx2fSums(aMatrixExtended, bVectorExtended, chi2sum,
geoIdVector.size(), scatteringMap, trackState,
logger);

geoIdVector.emplace_back(geoId);
}
}
}

/// @brief Calculate and update the covariance of the fitted parameters
///
/// This function calculates the covariance of the fitted parameters using
Expand Down Expand Up @@ -1185,85 +1283,15 @@ class Gx2Fitter {
Eigen::VectorXd bVectorExtended =
Eigen::VectorXd::Zero(dimsExtendedParams);

std::vector<BoundMatrix> jacobianFromStart;
jacobianFromStart.emplace_back(BoundMatrix::Identity());

// This vector stores the IDs for each visited material in order. We use
// it later for updating the scattering angles. We cannot use
// scatteringMap directly, since we cannot guarantee, that we will visit
// all stored material in each propagation.
std::vector<GeometryIdentifier> geoIdVector;

for (const auto& trackState : track.trackStates()) {
// Get and store geoId for the current surface
const GeometryIdentifier geoId =
trackState.referenceSurface().geometryId();
ACTS_DEBUG("Start to investigate trackState on surface " << geoId);
const auto typeFlags = trackState.typeFlags();
const bool stateHasMeasurement =
typeFlags.test(TrackStateFlag::MeasurementFlag);
const bool stateHasMaterial =
typeFlags.test(TrackStateFlag::MaterialFlag);

// First we figure out, if we would need to look into material surfaces
// at all. Later, we also check, if the material slab is valid,
// otherwise we modify this flag to ignore the material completely.
bool doMaterial = multipleScattering && stateHasMaterial;
if (doMaterial) {
const auto scatteringMapId = scatteringMap.find(geoId);
assert(scatteringMapId != scatteringMap.end() &&
"No scattering angles found for material surface.");
doMaterial = doMaterial && scatteringMapId->second.materialIsValid();
}

// We only consider states with a measurement (and/or material)
if (!stateHasMeasurement && !doMaterial) {
ACTS_DEBUG(" Skip state.");
continue;
}

// update all Jacobians from start
for (auto& jac : jacobianFromStart) {
jac = trackState.jacobian() * jac;
}

// Handle measurement
if (stateHasMeasurement) {
ACTS_DEBUG(" Handle measurement.");

const auto measDim = trackState.calibratedSize();

if (measDim < 1 || 6 < measDim) {
ACTS_ERROR("Can not process state with measurement with "
<< measDim << " dimensions.");
throw std::domain_error(
"Found measurement with less than 1 or more than 6 "
"dimension(s).");
}

countNdf += measDim;

visit_measurement(measDim, [&](auto N) {
addMeasurementToGx2fSums<N>(aMatrixExtended, bVectorExtended,
chi2sum, jacobianFromStart, trackState,
*m_addToSumLogger);
});
}

// Handle material
if (doMaterial) {
ACTS_DEBUG(" Handle material");
// Add for this material a new Jacobian, starting from this surface.
jacobianFromStart.emplace_back(BoundMatrix::Identity());

// Add the material contribution to the system
addMaterialToGx2fSums(aMatrixExtended, bVectorExtended, chi2sum,
geoIdVector.size(), scatteringMap, trackState,
*m_addToSumLogger);

geoIdVector.emplace_back(geoId);
}
}
fillGx2fSystem(track, aMatrixExtended, bVectorExtended, chi2sum, countNdf,
multipleScattering, scatteringMap, geoIdVector,
*m_addToSumLogger);

// Get required number of degrees of freedom ndfSystem.
// We have only 3 cases, because we always have l0, l1, phi, theta
Expand Down
Loading