StraightLineStepper.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file StraightLineStepper.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once

// Workaround for building on clang+libstdc++
#include "Acts/Utilities/detail/ReferenceWrapperAnyCompat.hpp"

#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Definitions/Direction.hpp"
#include "Acts/Definitions/Tolerance.hpp"
#include "Acts/Definitions/TrackParametrization.hpp"
#include "Acts/EventData/TrackParameters.hpp"
#include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/MagneticField/MagneticFieldContext.hpp"
#include "Acts/MagneticField/NullBField.hpp"
#include "Acts/Propagator/ConstrainedStep.hpp"
#include "Acts/Propagator/detail/SteppingHelper.hpp"
#include "Acts/Surfaces/BoundaryCheck.hpp"
#include "Acts/Surfaces/Surface.hpp"
#include "Acts/Utilities/Intersection.hpp"
#include "Acts/Utilities/Logger.hpp"
#include "Acts/Utilities/Result.hpp"

#include <algorithm>
#include <cmath>
#include <functional>
#include <limits>
#include <string>
#include <tuple>

namespace Acts {

class StraightLineStepper {
 public:
  using Jacobian = BoundMatrix;
  using Covariance = BoundSquareMatrix;
  using BoundState = std::tuple<BoundTrackParameters, Jacobian, double>;
  using CurvilinearState =
      std::tuple<CurvilinearTrackParameters, Jacobian, double>;
  using BField = NullBField;

  struct State {
    State() = delete;

    explicit State(const GeometryContext& gctx,
                   const MagneticFieldContext& mctx,
                   const BoundTrackParameters& par,
                   double ssize = std::numeric_limits<double>::max(),
                   double stolerance = s_onSurfaceTolerance)
        : particleHypothesis(par.particleHypothesis()),
          stepSize(ssize),
          tolerance(stolerance),
          geoContext(gctx) {
      (void)mctx;
      pars.template segment<3>(eFreePos0) = par.position(gctx);
      pars.template segment<3>(eFreeDir0) = par.direction();
      pars[eFreeTime] = par.time();
      pars[eFreeQOverP] = par.parameters()[eBoundQOverP];
      if (par.covariance()) {
        // Get the reference surface for navigation
        const auto& surface = par.referenceSurface();
        // set the covariance transport flag to true and copy
        covTransport = true;
        cov = BoundSquareMatrix(*par.covariance());
        jacToGlobal = surface.boundToFreeJacobian(gctx, par.parameters());
      }
    }

    BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();

    FreeMatrix jacTransport = FreeMatrix::Identity();

    Jacobian jacobian = Jacobian::Identity();

    FreeVector derivative = FreeVector::Zero();

    FreeVector pars = FreeVector::Zero();

    ParticleHypothesis particleHypothesis = ParticleHypothesis::pion();

    bool covTransport = false;
    Covariance cov = Covariance::Zero();

    double pathAccumulated = 0.;

    ConstrainedStep stepSize;

    // Previous step size for overstep estimation (ignored for SL stepper)
    double previousStepSize = 0.;

    double tolerance = s_onSurfaceTolerance;

    // Cache the geometry context of this propagation
    std::reference_wrapper<const GeometryContext> geoContext;
  };

  using state_type = State;

  StraightLineStepper() = default;

  State makeState(std::reference_wrapper<const GeometryContext> gctx,
                  std::reference_wrapper<const MagneticFieldContext> mctx,
                  const BoundTrackParameters& par,
                  double ssize = std::numeric_limits<double>::max(),
                  double stolerance = s_onSurfaceTolerance) const {
    return State{gctx, mctx, par, ssize, stolerance};
  }

  void resetState(
      State& state, const BoundVector& boundParams,
      const BoundSquareMatrix& cov, const Surface& surface,
      const double stepSize = std::numeric_limits<double>::max()) const;

  Result<Vector3> getField(State& state, const Vector3& pos) const {
    (void)state;
    (void)pos;
    // get the field from the cell
    return Result<Vector3>::success({0., 0., 0.});
  }

  Vector3 position(const State& state) const {
    return state.pars.template segment<3>(eFreePos0);
  }

  Vector3 direction(const State& state) const {
    return state.pars.template segment<3>(eFreeDir0);
  }

  double qOverP(const State& state) const { return state.pars[eFreeQOverP]; }

  double absoluteMomentum(const State& state) const {
    return particleHypothesis(state).extractMomentum(qOverP(state));
  }

  Vector3 momentum(const State& state) const {
    return absoluteMomentum(state) * direction(state);
  }

  double charge(const State& state) const {
    return particleHypothesis(state).extractCharge(qOverP(state));
  }

  const ParticleHypothesis& particleHypothesis(const State& state) const {
    return state.particleHypothesis;
  }

  double time(const State& state) const { return state.pars[eFreeTime]; }

  double overstepLimit(const State& state) const {
    (void)state;
    return -m_overstepLimit;
  }

  Intersection3D::Status updateSurfaceStatus(
      State& state, const Surface& surface, Direction navDir,
      const BoundaryCheck& bcheck,
      ActsScalar surfaceTolerance = s_onSurfaceTolerance,
      const Logger& logger = getDummyLogger()) const {
    return detail::updateSingleSurfaceStatus<StraightLineStepper>(
        *this, state, surface, navDir, bcheck, surfaceTolerance, logger);
  }

  template <typename object_intersection_t>
  void updateStepSize(State& state, const object_intersection_t& oIntersection,
                      Direction /*direction*/, bool release = true) const {
    detail::updateSingleStepSize<StraightLineStepper>(state, oIntersection,
                                                      release);
  }

  void setStepSize(State& state, double stepSize,
                   ConstrainedStep::Type stype = ConstrainedStep::actor,
                   bool release = true) const {
    state.previousStepSize = state.stepSize.value();
    state.stepSize.update(stepSize, stype, release);
  }

  double getStepSize(const State& state, ConstrainedStep::Type stype) const {
    return state.stepSize.value(stype);
  }

  void releaseStepSize(State& state) const {
    state.stepSize.release(ConstrainedStep::actor);
  }

  std::string outputStepSize(const State& state) const {
    return state.stepSize.toString();
  }

  Result<BoundState> boundState(
      State& state, const Surface& surface, bool transportCov = true,
      const FreeToBoundCorrection& freeToBoundCorrection =
          FreeToBoundCorrection(false)) const;

  CurvilinearState curvilinearState(State& state,
                                    bool transportCov = true) const;

  void update(State& state, const FreeVector& freeParams,
              const BoundVector& boundParams, const Covariance& covariance,
              const Surface& surface) const;

  void update(State& state, const Vector3& uposition, const Vector3& udirection,
              double qop, double time) const;

  void transportCovarianceToCurvilinear(State& state) const;

  void transportCovarianceToBound(
      State& state, const Surface& surface,
      const FreeToBoundCorrection& freeToBoundCorrection =
          FreeToBoundCorrection(false)) const;

  template <typename propagator_state_t, typename navigator_t>
  Result<double> step(propagator_state_t& state,
                      const navigator_t& /*navigator*/) const {
    // use the adjusted step size
    const auto h = state.stepping.stepSize.value() * state.options.direction;
    const auto m = state.stepping.particleHypothesis.mass();
    const auto p = absoluteMomentum(state.stepping);
    // time propagates along distance as 1/b = sqrt(1 + m²/p²)
    const auto dtds = std::hypot(1., m / p);
    // Update the track parameters according to the equations of motion
    Vector3 dir = direction(state.stepping);
    state.stepping.pars.template segment<3>(eFreePos0) += h * dir;
    state.stepping.pars[eFreeTime] += h * dtds;
    // Propagate the jacobian
    if (state.stepping.covTransport) {
      // The step transport matrix in global coordinates
      FreeMatrix D = FreeMatrix::Identity();
      D.block<3, 3>(0, 4) = ActsSquareMatrix<3>::Identity() * h;
      // Extend the calculation by the time propagation
      // Evaluate dt/dlambda
      D(3, 7) = h * m * m * state.stepping.pars[eFreeQOverP] / dtds;
      // Set the derivative factor the time
      state.stepping.derivative(3) = dtds;
      // Update jacobian and derivative
      state.stepping.jacTransport = D * state.stepping.jacTransport;
      state.stepping.derivative.template head<3>() = dir;
    }
    // state the path length
    state.stepping.pathAccumulated += h;

    // return h
    return h;
  }

 private:
  double m_overstepLimit = s_onSurfaceTolerance;
};

}  // namespace Acts