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/EventData/TrackParameters.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/Surface.hpp"
#include "Acts/Utilities/Definitions.hpp"
#include "Acts/Utilities/Intersection.hpp"
#include "Acts/Utilities/Result.hpp"

#include <cmath>
#include <functional>

namespace Acts {

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

  struct State {
    State() = delete;

    template <typename parameters_t>
    explicit State(std::reference_wrapper<const GeometryContext> gctx,
                   std::reference_wrapper<const MagneticFieldContext> /*mctx*/,
                   const parameters_t& par, NavigationDirection ndir = forward,
                   double ssize = std::numeric_limits<double>::max(),
                   double stolerance = s_onSurfaceTolerance)
        : pos(par.position(gctx)),
          dir(par.unitDirection()),
          p(par.absoluteMomentum()),
          q(par.charge()),
          t(par.time()),
          navDir(ndir),
          stepSize(ndir * std::abs(ssize)),
          tolerance(stolerance),
          geoContext(gctx) {
      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 = BoundSymMatrix(*par.covariance());
        surface.initJacobianToGlobal(gctx, jacToGlobal, pos, dir,
                                     par.parameters());
      }
    }

    BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();

    FreeMatrix jacTransport = FreeMatrix::Identity();

    Jacobian jacobian = Jacobian::Identity();

    FreeVector derivative = FreeVector::Zero();

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

    Vector3D pos = Vector3D(0., 0., 0.);

    Vector3D dir = Vector3D(1., 0., 0.);

    double p = 0.;

    double q = 0.;

    double t = 0.;

    NavigationDirection navDir;

    double pathAccumulated = 0.;

    ConstrainedStep stepSize = std::numeric_limits<double>::max();

    // 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;

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

  Vector3D getField(State& /*state*/, const Vector3D& /*pos*/) const {
    // get the field from the cell
    return Vector3D(0., 0., 0.);
  }

  Vector3D position(const State& state) const { return state.pos; }

  Vector3D direction(const State& state) const { return state.dir; }

  double momentum(const State& state) const { return state.p; }

  double charge(const State& state) const { return state.q; }

  double time(const State& state) const { return state.t; }

  double overstepLimit(const State& /*state*/) const {
    return s_onSurfaceTolerance;
  }

  Intersection3D::Status updateSurfaceStatus(
      State& state, const Surface& surface, const BoundaryCheck& bcheck) const {
    return detail::updateSingleSurfaceStatus<StraightLineStepper>(
        *this, state, surface, bcheck);
  }

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

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

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

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

  BoundState boundState(State& state, const Surface& surface) const;

  CurvilinearState curvilinearState(State& state) const;

  void update(State& state, const FreeVector& parameters,
              const Covariance& covariance) const;

  void update(State& state, const Vector3D& uposition,
              const Vector3D& udirection, double up, double time) const;

  void covarianceTransport(State& state) const;

  void covarianceTransport(State& state, const Surface& surface) const;

  template <typename propagator_state_t>
  Result<double> step(propagator_state_t& state) const {
    // use the adjusted step size
    const auto h = state.stepping.stepSize;
    // time propagates along distance as 1/b = sqrt(1 + m²/p²)
    const auto dtds = std::hypot(1., state.options.mass / state.stepping.p);
    // Update the track parameters according to the equations of motion
    state.stepping.pos += h * state.stepping.dir;
    state.stepping.t += 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) = ActsSymMatrixD<3>::Identity() * h;
      // Extend the calculation by the time propagation
      // Evaluate dt/dlambda
      D(3, 7) = h * state.options.mass * state.options.mass *
                (state.stepping.q == 0. ? 1. : state.stepping.q) /
                (state.stepping.p * 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>() = state.stepping.dir;
    }
    // state the path length
    state.stepping.pathAccumulated += h;

    // return h
    return h;
  }
};

}  // namespace Acts