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-2022 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/Tolerance.hpp"
#include "Acts/Definitions/Units.hpp"
#include "Acts/EventData/TrackParameters.hpp"
#include "Acts/MagneticField/MagneticFieldProvider.hpp"
#include "Acts/Propagator/ConstrainedStep.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 <cmath>
#include <functional>
#include <limits>

namespace Acts {

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

  struct State {
    State() = delete;

    explicit State(const GeometryContext& gctx,
                   MagneticFieldProvider::Cache fieldCacheIn,
                   const BoundTrackParameters& par,
                   double ssize = std::numeric_limits<double>::max())
        : particleHypothesis(par.particleHypothesis()),
          stepSize(ssize),
          fieldCache(std::move(fieldCacheIn)),
          geoContext(gctx) {
      pars.template segment<3>(eFreePos0) = par.position(gctx);
      pars.template segment<3>(eFreeDir0) = par.direction();
      pars[eFreeTime] = par.time();
      pars[eFreeQOverP] = par.parameters()[eBoundQOverP];

      // 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 = BoundSquareMatrix(*par.covariance());
        jacToGlobal = surface.boundToFreeJacobian(gctx, par.parameters());
      }
    }

    FreeVector pars = FreeVector::Zero();

    ParticleHypothesis particleHypothesis = ParticleHypothesis::pion();

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

    Jacobian jacobian = Jacobian::Identity();

    BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();

    FreeMatrix jacTransport = FreeMatrix::Identity();

    FreeVector derivative = FreeVector::Zero();

    double pathAccumulated = 0.;

    ConstrainedStep stepSize;

    double previousStepSize = 0.;

    MagneticFieldProvider::Cache fieldCache;

    std::reference_wrapper<const GeometryContext> geoContext;

    extensionlist_t extension;

    auctioneer_t auctioneer;

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

  EigenStepper(std::shared_ptr<const MagneticFieldProvider> bField,
               double overstepLimit = 100 * UnitConstants::um);

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

  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 {
    // get the field from the cell
    return m_bField->getField(pos, state.fieldCache);
  }

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

  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<EigenStepper>(
        *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<EigenStepper>(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();
  }

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

  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 qOverP, 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;

  void setIdentityJacobian(State& state) const;

 protected:
  std::shared_ptr<const MagneticFieldProvider> m_bField;

  double m_overstepLimit;
};
}  // namespace Acts

#include "Acts/Propagator/EigenStepper.ipp"