Propagator.hpp

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

// This file is part of the Acts project.
//
// Copyright (C) 2016-2019 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 <cmath>
#include <functional>
#include <memory>
#include <type_traits>

#include <boost/algorithm/string.hpp>

#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/MagneticField/MagneticFieldContext.hpp"
#include "Acts/Propagator/AbortList.hpp"
#include "Acts/Propagator/ActionList.hpp"
#include "Acts/Propagator/PropagatorError.hpp"
#include "Acts/Propagator/StandardAborters.hpp"
#include "Acts/Propagator/StepperConcept.hpp"
#include "Acts/Propagator/detail/LoopProtection.hpp"
#include "Acts/Propagator/detail/VoidPropagatorComponents.hpp"
#include "Acts/Utilities/Definitions.hpp"
#include "Acts/Utilities/Result.hpp"
#include "Acts/Utilities/Units.hpp"

namespace Acts {

template <typename parameters_t, typename... result_list>
struct PropagatorResult : private detail::Extendable<result_list...> {
  using detail::Extendable<result_list...>::get;

  std::unique_ptr<const parameters_t> endParameters = nullptr;

  std::unique_ptr<const BoundMatrix> transportJacobian = nullptr;

  unsigned int steps = 0;

  double pathLength = 0.;
};

template <typename action_list_t = ActionList<>,
          typename aborter_list_t = AbortList<>>
struct PropagatorOptions {
  using action_list_type = action_list_t;
  using aborter_list_type = aborter_list_t;

  PropagatorOptions() = delete;

  PropagatorOptions(
      const PropagatorOptions<action_list_t, aborter_list_t>& po) = default;

  PropagatorOptions(std::reference_wrapper<const GeometryContext> gctx,
                    std::reference_wrapper<const MagneticFieldContext> mctx)
      : geoContext(gctx), magFieldContext(mctx) {}

  template <typename extended_aborter_list_t>
  PropagatorOptions<action_list_t, extended_aborter_list_t> extend(
      extended_aborter_list_t aborters) const {
    PropagatorOptions<action_list_t, extended_aborter_list_t> eoptions(
        geoContext, magFieldContext);
    // Copy the options over
    eoptions.direction = direction;
    eoptions.absPdgCode = absPdgCode;
    eoptions.mass = mass;
    eoptions.maxSteps = maxSteps;
    eoptions.maxRungeKuttaStepTrials = maxRungeKuttaStepTrials;
    eoptions.maxStepSize = direction * std::abs(maxStepSize);
    eoptions.targetTolerance = targetTolerance;
    eoptions.pathLimit = direction * std::abs(pathLimit);
    eoptions.loopProtection = loopProtection;
    eoptions.loopFraction = loopFraction;
    // Output option
    eoptions.debug = debug;
    eoptions.debugString = debugString;
    eoptions.debugPfxWidth = debugPfxWidth;
    eoptions.debugMsgWidth = debugMsgWidth;
    // Stepper options
    eoptions.tolerance = tolerance;
    eoptions.stepSizeCutOff = stepSizeCutOff;
    // Action / abort list
    eoptions.actionList = std::move(actionList);
    eoptions.abortList = std::move(aborters);
    // And return the options
    return eoptions;
  }

  NavigationDirection direction = forward;

  int absPdgCode = 211;

  double mass = 139.57018 * UnitConstants::MeV;

  unsigned int maxSteps = 1000;

  unsigned int maxRungeKuttaStepTrials = 10000;

  double maxStepSize = std::numeric_limits<double>::max();

  double pathLimit = std::numeric_limits<double>::max();

  double targetTolerance = s_onSurfaceTolerance;

  bool loopProtection = true;
  double loopFraction = 0.5;  

  //  -> @todo: move to a debug struct
  // - the string where debug messages are stored (optionally)
  // - it also has some formatting options
  bool debug = false;            
  std::string debugString = "";  
  size_t debugPfxWidth = 30;     
  size_t debugMsgWidth = 50;     

  // Configurations for Stepper
  double tolerance = 1e-4;

  double stepSizeCutOff = 0.;

  action_list_t actionList;

  aborter_list_t abortList;

  std::reference_wrapper<const GeometryContext> geoContext;

  std::reference_wrapper<const MagneticFieldContext> magFieldContext;
};

template <typename stepper_t, typename navigator_t = detail::VoidNavigator>
class Propagator final {
  using Jacobian = BoundMatrix;
  using BoundState = std::tuple<BoundParameters, Jacobian, double>;
  using CurvilinearState = std::tuple<CurvilinearParameters, Jacobian, double>;

  static_assert(StepperStateConcept<typename stepper_t::State>,
                "Stepper does not fulfill stepper concept.");
  static_assert(StepperConcept<stepper_t>,
                "Stepper does not fulfill stepper concept.");

 public:
  using Stepper = stepper_t;

  using Navigator = navigator_t;

  using StepperState = typename Stepper::State;

  using NavigatorState = typename navigator_t::State;

  explicit Propagator(stepper_t stepper, navigator_t navigator = navigator_t())
      : m_stepper(std::move(stepper)), m_navigator(std::move(navigator)) {}

  template <typename propagator_options_t>
  struct State {
    template <typename parameters_t>
    State(const parameters_t& start, const propagator_options_t& topts)
        : options(topts),
          stepping(topts.geoContext, topts.magFieldContext, start,
                   topts.direction, topts.maxStepSize, topts.tolerance),
          geoContext(topts.geoContext) {
      // Setting the start surface
      navigation.startSurface = &start.referenceSurface();
    }

    propagator_options_t options;

    StepperState stepping;

    NavigatorState navigation;

    std::reference_wrapper<const GeometryContext> geoContext;
  };

 private:
  template <typename parameters_t, typename action_list_t>
  struct result_type_helper {
    template <typename... args>
    using this_result_type = PropagatorResult<parameters_t, args...>;

    using type = typename action_list_t::template result_type<this_result_type>;
  };

  template <typename parameters_t, typename action_list_t>
  using action_list_t_result_t =
      typename result_type_helper<parameters_t, action_list_t>::type;

  template <typename result_t, typename propagator_state_t>
  Result<result_t> propagate_impl(propagator_state_t& state) const;

 public:
  template <typename parameters_t, typename propagator_options_t,
            typename path_aborter_t = PathLimitReached>
  Result<action_list_t_result_t<
      CurvilinearParameters, typename propagator_options_t::action_list_type>>
  propagate(const parameters_t& start,
            const propagator_options_t& options) const;

  template <typename parameters_t, typename propagator_options_t,
            typename target_aborter_t = SurfaceReached,
            typename path_aborter_t = PathLimitReached>
  Result<action_list_t_result_t<
      BoundParameters, typename propagator_options_t::action_list_type>>
  propagate(const parameters_t& start, const Surface& target,
            const propagator_options_t& options) const;

 private:
  stepper_t m_stepper;

  navigator_t m_navigator;

  template <typename propagator_state_t>
  void debugLog(propagator_state_t& state,
                const std::function<std::string()>& logAction) const;
};

}  // namespace Acts

#include "Acts/Propagator/Propagator.ipp"