Track_fitting

Namespaces
namespace	MuonExperimental::Gx2fConstants

Classes
struct	MuonExperimental::Gx2fMaterialProperties
	A container to store material properties for each material surface. More...
struct	MuonExperimental::Gx2FitterExtensions< traj_t >
	Extension struct which holds delegates to customise the GX2F behaviour. More...
struct	MuonExperimental::Gx2FitterOptions< traj_t >
	Combined options for the Global-Chi-Square fitter. More...
struct	MuonExperimental::Gx2FitterResult< traj_t >
	Result container for a global chi-square fit. More...
struct	MuonExperimental::Gx2fSystem
	A container to manage all properties of a gx2f system. More...
class	MuonExperimental::Gx2Fitter< propagator_t, traj_t >
	Global Chi Square fitter (GX2F) implementation. More...

Functions
void	MuonExperimental::addMeasurementToGx2fSumsBackend (Gx2fSystem &extendedSystem, const std::vector< BoundMatrix > &jacobianFromStart, const std::vector< std::size_t > &materialIndices, const std::vector< double > &qOverPderiv, const Eigen::MatrixXd &covarianceMeasurement, const BoundVector &predicted, const Eigen::VectorXd &measurement, const Eigen::MatrixXd &projector, const Logger &logger)
	Adds a measurement to the GX2F equation system in a modular backend function.
template<std::size_t kMeasDim, typename track_state_t>
void	MuonExperimental::addMeasurementToGx2fSums (Gx2fSystem &extendedSystem, const std::vector< BoundMatrix > &jacobianFromStart, const std::vector< std::size_t > &materialIndices, const std::vector< double > &qOverPderiv, const track_state_t &trackState, const Logger &logger)
	Process measurements and fill the aMatrix and bVector.
template<TrackProxyConcept track_proxy_t>
void	MuonExperimental::fillGx2fSystem (const track_proxy_t track, Gx2fSystem &extendedSystem, const std::unordered_map< GeometryIdentifier, Gx2fMaterialProperties > &materialMap, const std::vector< std::size_t > &materialParamIndices, std::vector< GeometryIdentifier > &geoIdVector, const Logger &logger)
	Fill the GX2F system with data from a track.
Eigen::VectorXd	MuonExperimental::computeGx2fDeltaParams (const Gx2fSystem &extendedSystem)
	Solve the gx2f system to get the delta parameters for the update.
void	MuonExperimental::updateGx2fParams (BoundTrackParameters &params, const Eigen::VectorXd &deltaParamsExtended, std::unordered_map< GeometryIdentifier, Gx2fMaterialProperties > &materialMap, const std::vector< GeometryIdentifier > &geoIdVector)
	Update parameters (and scattering angles if applicable).
void	MuonExperimental::updateGx2fCovarianceParams (BoundMatrix &fullCovariancePredicted, Gx2fSystem &extendedSystem)
	Calculate and update the covariance of the fitted parameters.

Detailed Description

Function Documentation

addMeasurementToGx2fSums()

template<std::size_t kMeasDim, typename track_state_t>

void MuonExperimental::addMeasurementToGx2fSums	(	Gx2fSystem &	extendedSystem,
		const std::vector< BoundMatrix > &	jacobianFromStart,
		const std::vector< std::size_t > &	materialIndices,
		const std::vector< double > &	qOverPderiv,
		const track_state_t &	trackState,
		const Logger &	logger )

Process measurements and fill the aMatrix and bVector.

The function processes each measurement for the GX2F Actor fitting process. It extracts the information from the track state and adds it to aMatrix, bVector, and chi2sum.

Template Parameters

kMeasDim	Number of dimensions of the measurement
track_state_t	The type of the track state

Parameters

extendedSystem	All parameters of the current equation system to update
jacobianFromStart	The Jacobian matrix from start to the current state
trackState	The track state to analyse
logger	A logger instance

Definition at line 545 of file Chi2Fitter.h.

                                                    {
  const SquareMatrix<kMeasDim> covarianceMeasurement =
      trackState.template calibratedCovariance<kMeasDim>();
 
  const BoundVector predicted = trackState.smoothed();
 
  const Vector<kMeasDim> measurement =
      trackState.template calibrated<kMeasDim>();
 
  const Matrix<kMeasDim, eBoundSize> projector =
      trackState.template projectorSubspaceHelper<kMeasDim>().projector();
 
  addMeasurementToGx2fSumsBackend(
      extendedSystem, jacobianFromStart, materialIndices, qOverPderiv,
      covarianceMeasurement, predicted, measurement, projector, logger);
}

addMeasurementToGx2fSumsBackend()

void MuonExperimental::addMeasurementToGx2fSumsBackend	(	Gx2fSystem &	extendedSystem,
		const std::vector< BoundMatrix > &	jacobianFromStart,
		const std::vector< std::size_t > &	materialIndices,
		const std::vector< double > &	qOverPderiv,
		const Eigen::MatrixXd &	covarianceMeasurement,
		const BoundVector &	predicted,
		const Eigen::VectorXd &	measurement,
		const Eigen::MatrixXd &	projector,
		const Logger &	logger )

Adds a measurement to the GX2F equation system in a modular backend function.

This function processes measurement data and integrates it into the GX2F system.

Parameters

extendedSystem	All parameters of the current equation system to update.
jacobianFromStart	The Jacobian matrix from the start to the current state.
covarianceMeasurement	The covariance matrix of the measurement.
predicted	The predicted state vector based on the track state.
measurement	The measurement vector.
projector	The projection matrix.
logger	A logger instance.

Note

The dynamic Eigen matrices are suboptimal. We could think of templating again in the future on kMeasDims. We currently use dynamic matrices to reduce the memory during compile time.

Definition at line 213 of file Chi2Fitter.cxx.

                          {
  // First, we try to invert the covariance matrix. If the inversion fails, we
  // can already abort.
  const auto safeInvCovMeasurement = safeInverse(covarianceMeasurement);
  if (!safeInvCovMeasurement) {
    ACTS_WARNING("addMeasurementToGx2fSums: safeInvCovMeasurement failed.");
    ACTS_VERBOSE("    covarianceMeasurement:\n" << covarianceMeasurement);
    return;
  }
  ACTS_VERBOSE("Jacobian size: " << jacobianFromStart.size());
  // Create an extended Jacobian. This one contains only eBoundSize rows,
  // because the rest is irrelevant. We fill it in the next steps.
  // TODO make dimsExtendedParams template with unrolling
  Eigen::MatrixXd extendedJacobian =
      Eigen::MatrixXd::Zero(eBoundSize, extendedSystem.nDims());
 
  // This part of the Jacobian comes from the material-less propagation
  extendedJacobian.topLeftCorner<eBoundSize, eBoundSize>() =
      jacobianFromStart[0];
 
  // If we have material, loop here over all Jacobians. We add extra columns for
  // their phi-theta projections. These parts account for the propagation of the
  // scattering angles. if energy loss is enabled, we also add the q/p
  // projection.
  for (std::size_t matSurface = 1; matSurface < jacobianFromStart.size();
       matSurface++) {
    const BoundMatrix& jac = jacobianFromStart[matSurface];
 
    // The position, where we need to insert the values in the extended Jacobian
    // accounting for material effects
    const std::size_t deltaPosition = materialIndices.at(matSurface - 1);
    ACTS_DEBUG("Update material jacobian " << deltaPosition << ", "
                                           << extendedSystem.nDims());
 
    // check the dimension of the material map: 2 means only scattering, 3 means
    // scattering and energy loss amd 1 only energy loss
    const std::size_t matDim = materialIndices.at(matSurface) - deltaPosition;
 
    // check if energy loss is enabled, if yes, we need to add the q/p
    // projection in the extended jacobian
    switch (matDim) {
      case 1: {  // energy loss only
        Acts::Matrix<eBoundSize, 1> proj = qOverPProjector;
        proj(eBoundQOverP, 0) = qOverPderiv.at(matSurface - 1);
        extendedJacobian.block<eBoundSize, 1>(0, deltaPosition) = jac * proj;
        break;
      }
      case 2:  // scattering only
        extendedJacobian.block<eBoundSize, 2>(0, deltaPosition) =
            jac * phiThetaProjector;
        break;
      case 3: {  // scattering and energy loss
        Acts::Matrix<eBoundSize, 3> proj = phiThetaQOverPProjector;
        proj(eBoundQOverP, 2) = qOverPderiv.at(matSurface - 1);
        extendedJacobian.block<eBoundSize, 3>(0, deltaPosition) = jac * proj;
        break;
      }
      default:
        ACTS_ERROR("Invalid material dimension: "
                   << matDim << " - this should be 1, 2, or 3.");
        throw std::domain_error("Invalid material dimension.");
    }
  }
 
  const Eigen::MatrixXd projJacobian = projector * extendedJacobian;
 
  const Eigen::VectorXd projPredicted = projector * predicted;
 
  const Eigen::VectorXd residual = measurement - projPredicted;
 
  // Finally contribute to chi2sum, aMatrix, and bVector
  extendedSystem.chi2() +=
      (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0);
 
  extendedSystem.aMatrix() +=
      (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian)
          .eval();
 
  extendedSystem.bVector() +=
      (residual.transpose() * (*safeInvCovMeasurement) * projJacobian)
          .eval()
          .transpose();
 
  ACTS_VERBOSE(
      "Contributions in addMeasurementToGx2fSums:\n"
      << "    predicted:   " << predicted.transpose() << "\n"
      << "    measurement: " << measurement.transpose() << "\n"
      << "    covarianceMeasurement:\n"
      << covarianceMeasurement << "\n"
      << "    projector:\n"
      << projector.eval() << "\n"
      << "    projJacobian:\n"
      << projJacobian.eval() << "\n"
      << "    projPredicted: " << (projPredicted.transpose()).eval() << "\n"
      << "    residual: " << (residual.transpose()).eval() << "\n"
      << "    extendedJacobian:\n"
      << extendedJacobian << "\n"
      << "    aMatrix contribution:\n"
      << (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian)
             .eval()
      << "\n"
      << "    bVector contribution: "
      << (residual.transpose() * (*safeInvCovMeasurement) * projJacobian).eval()
      << "\n"
      << "    chi2sum contribution: "
      << (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0)
      << "\n"
      << "    safeInvCovMeasurement:\n"
      << (*safeInvCovMeasurement));
}

computeGx2fDeltaParams()

Eigen::VectorXd MuonExperimental::computeGx2fDeltaParams

(

const Gx2fSystem &

extendedSystem

)

Solve the gx2f system to get the delta parameters for the update.

This function computes the delta parameters for the GX2F Actor fitting process by solving the linear equation system [a] * delta = b. It uses the column-pivoting Householder QR decomposition for numerical stability.

Parameters

extendedSystem

All parameters of the current equation system

Returns

Delta parameters for the GX2F update

Definition at line 331 of file Chi2Fitter.cxx.

                                                                       {
  return extendedSystem.aMatrix().colPivHouseholderQr().solve(
      extendedSystem.bVector());
}

fillGx2fSystem()

template<TrackProxyConcept track_proxy_t>

void MuonExperimental::fillGx2fSystem	(	const track_proxy_t	track,
		Gx2fSystem &	extendedSystem,
		const std::unordered_map< GeometryIdentifier, Gx2fMaterialProperties > &	materialMap,
		const std::vector< std::size_t > &	materialParamIndices,
		std::vector< GeometryIdentifier > &	geoIdVector,
		const Logger &	logger )

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.

Template Parameters

track_proxy_t

The type of the track proxy

Parameters

track	A constant track proxy to inspect
extendedSystem	All parameters of the current equation system
multipleScattering	Flag to consider multiple scattering in the calculation
materialMap	Map of geometry identifiers to scattering properties, containing scattering angles and validation status
logger	A logger instance

Invalid material state

Definition at line 582 of file Chi2Fitter.h.

                                                                      {
  std::vector<BoundMatrix> jacobianFromStart;
  jacobianFromStart.emplace_back(BoundMatrix::Identity());
  geoIdVector.reserve(track.nTrackStates());
  std::vector<double> qOverPderiv;
 
  bool foundELoss{false};
  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();
 
    // update all Jacobians from start for every track state
    for (auto& jac : jacobianFromStart) {
      jac = trackState.jacobian() * jac;
    }
 
    // Skip the hole states
    if (!typeFlags.hasMeasurement() && !typeFlags.hasMaterial()) {
      ACTS_DEBUG("    Skip state.");
      continue;
    }
 
    // Handle measurement
    if (typeFlags.hasMeasurement()) {
      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).");
      }
 
      extendedSystem.ndf() += measDim;
 
      visit_measurement(measDim, [&](auto N) {
        addMeasurementToGx2fSums<N>(extendedSystem, jacobianFromStart,
                                    materialParamIndices, qOverPderiv,
                                    trackState, logger);
      });
    }
 
    // Handle material
    if (typeFlags.hasMaterial()) {
      ACTS_DEBUG("    Handle material");
      const auto materialMapId = materialMap.find(geoId);
      if (materialMapId == materialMap.end()) {
        ACTS_ERROR("No material properties found for surface " << geoId);
        throw std::runtime_error("No material properties found for surface.");
      }
      if (materialMapId->second.nDim() == 0ul) {
        continue;
      }
      if (materialMapId->second.energyLoss().isValid()) {
        ACTS_VERBOSE("The material surface carries energy loss");
        foundELoss = true;
      }
      const std::size_t deltaPosition =
          materialParamIndices.at(geoIdVector.size());
      geoIdVector.push_back(geoId);
 
      materialMapId->second.contributionToGx2fSums(
          extendedSystem, trackState.smoothed()[eBoundTheta], deltaPosition,
          logger);
 
      // after material surfaces, the measurements depend on the updated qoverp
      // we need to scale the jacobian term for the qoverp row
      if (foundELoss) {
        ACTS_DEBUG("    Update Jacobian for q/p after material surface "
                   << geoId << " with ratio " << materialMapId->second.qOverPratio());
        jacobianFromStart.front().row(eBoundQOverP) *=
            materialMapId->second.qOverPratio();
      }
      // Add for this material a new Jacobian, starting from this surface.
      jacobianFromStart.emplace_back(BoundMatrix::Identity());
      qOverPderiv.push_back(materialMapId->second.qOverPderiv());
    }
  }
}

updateGx2fCovarianceParams()

void MuonExperimental::updateGx2fCovarianceParams	(	BoundMatrix &	fullCovariancePredicted,
		Gx2fSystem &	extendedSystem )

Calculate and update the covariance of the fitted parameters.

This function calculates the covariance of the fitted parameters using cov = inv([a]) It then updates the first square block of size ndfSystem. This ensures, that we only update the covariance for fitted parameters. (In case of no qop/time fit)

Parameters

fullCovariancePredicted	The covariance matrix to update
extendedSystem	All parameters of the current equation system

Definition at line 195 of file Chi2Fitter.cxx.

                                                            {
  const std::size_t nDof = extendedSystem.findRequiredNdf();
  // make invertible
  for (std::size_t i = 0; i < extendedSystem.nDims(); ++i) {
    if (extendedSystem.aMatrix()(i, i) == 0.) {
      extendedSystem.aMatrix()(i, i) = 1.;
    }
  }
 
  visit_measurement(nDof, [&](auto N) {
    fullCovariancePredicted.topLeftCorner<N, N>() =
        extendedSystem.aMatrix().inverse().topLeftCorner<N, N>();
  });
 
  return;
}

updateGx2fParams()

void MuonExperimental::updateGx2fParams	(	BoundTrackParameters &	params,
		const Eigen::VectorXd &	deltaParamsExtended,
		std::unordered_map< GeometryIdentifier, Gx2fMaterialProperties > &	materialMap,
		const std::vector< GeometryIdentifier > &	geoIdVector )

Update parameters (and scattering angles if applicable).

Parameters

params	Parameters to be updated
deltaParamsExtended	Delta parameters for bound parameter and scattering angles
materialMap	Map of geometry identifiers to scattering properties, containing all scattering angles and covariances
geoIdVector	Vector of geometry identifiers corresponding to material surfaces

Definition at line 175 of file Chi2Fitter.cxx.

                                                      {
  // update params
  params.parameters() +=
      deltaParamsExtended.topLeftCorner<eBoundSize, 1>().eval();
 
  std::size_t deltaPosition{eBoundSize};
  // update the scattering angles and energy loss.
  for (const GeometryIdentifier& geoId : geoIdVector) {
    const auto materialMapId = materialMap.find(geoId);
    assert(materialMapId != materialMap.end() &&
           "No material properties found for material surface.");
 
    materialMapId->second.updateParameters(deltaParamsExtended, deltaPosition);
    deltaPosition += materialMapId->second.nDim();
  }
}