EigenStepper.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EigenStepper.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/Propagator/DefaultExtension.hpp"
#include "Acts/Propagator/DenseEnvironmentExtension.hpp"
#include "Acts/Propagator/EigenStepperError.hpp"
#include "Acts/Propagator/StepperExtensionList.hpp"
#include "Acts/Propagator/detail/Auctioneer.hpp"
#include "Acts/Propagator/detail/SteppingHelper.hpp"
#include "Acts/Utilities/Intersection.hpp"
#include "Acts/Utilities/Result.hpp"
#include "Acts/Utilities/Units.hpp"

#include <cmath>
#include <functional>
#include <limits>

namespace Acts {

using namespace Acts::UnitLiterals;

template <typename bfield_t,
          typename extensionlist_t = StepperExtensionList<DefaultExtension>,
          typename auctioneer_t = detail::VoidAuctioneer>
class EigenStepper {
 public:
  using Jacobian = BoundMatrix;
  using Covariance = BoundSymMatrix;
  using BoundState = std::tuple<BoundTrackParameters, Jacobian, double>;
  using CurvilinearState =
      std::tuple<CurvilinearTrackParameters, Jacobian, double>;
  using BField = bfield_t;

  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),
          fieldCache(mctx),
          geoContext(gctx) {
      // Init the jacobian matrix if needed
      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());
      }
    }

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

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

    double p = 0.;

    double q = 1.;

    double t = 0.;

    NavigationDirection navDir;

    Jacobian jacobian = Jacobian::Identity();

    BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();

    FreeMatrix jacTransport = FreeMatrix::Identity();

    FreeVector derivative = FreeVector::Zero();

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

    double pathAccumulated = 0.;

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

    double previousStepSize = 0.;

    double tolerance = s_onSurfaceTolerance;

    typename BField::Cache fieldCache;

    std::reference_wrapper<const GeometryContext> geoContext;

    extensionlist_t extension;

    auctioneer_t auctioneer;

    struct {
      Vector3D B_first, B_middle, B_last;
      Vector3D k1, k2, k3, k4;
      std::array<double, 4> kQoP;
    } stepData;
  };

  EigenStepper(BField bField);

  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 m_bField.getField(pos, state.fieldCache);
  }

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

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

  template <typename object_intersection_t>
  void updateStepSize(State& state, const object_intersection_t& oIntersection,
                      bool release = true) const {
    detail::updateSingleStepSize<EigenStepper>(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();
  }

  double overstepLimit(const State& /*state*/) const {
    // A dynamic overstep limit could sit here
    return -m_overstepLimit;
  }

  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;

 private:
  BField m_bField;

  double m_overstepLimit = 100_um;
};
}  // namespace Acts

#include "Acts/Propagator/EigenStepper.ipp"