DefaultExtension.hpp

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

// This file is part of the Acts project.
//
// Copyright (C) 2018-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

#include <array>

#include "Acts/Utilities/Helpers.hpp"

namespace Acts {

struct DefaultExtension {
  DefaultExtension() = default;

  template <typename propagator_state_t, typename stepper_t>
  int bid(const propagator_state_t& /*unused*/,
          const stepper_t& /*unused*/) const {
    return 1;
  }

  template <typename propagator_state_t, typename stepper_t>
  bool k(const propagator_state_t& state, const stepper_t& stepper,
         Vector3D& knew, const Vector3D& bField, std::array<double, 4>& kQoP,
         const int i = 0, const double h = 0.,
         const Vector3D& kprev = Vector3D()) {
    auto qop =
        stepper.charge(state.stepping) / stepper.momentum(state.stepping);
    // First step does not rely on previous data
    if (i == 0) {
      knew = qop * stepper.direction(state.stepping).cross(bField);
      kQoP = {0., 0., 0., 0.};
    } else {
      knew =
          qop * (stepper.direction(state.stepping) + h * kprev).cross(bField);
    }
    return true;
  }

  template <typename propagator_state_t, typename stepper_t>
  bool finalize(propagator_state_t& state, const stepper_t& stepper,
                const double h) const {
    propagateTime(state, stepper, h);
    return true;
  }

  template <typename propagator_state_t, typename stepper_t>
  bool finalize(propagator_state_t& state, const stepper_t& stepper,
                const double h, FreeMatrix& D) const {
    propagateTime(state, stepper, h);
    return transportMatrix(state, stepper, h, D);
  }

 private:
  template <typename propagator_state_t, typename stepper_t>
  void propagateTime(propagator_state_t& state, const stepper_t& stepper,
                     const double h) const {
    auto derivative =
        std::hypot(1, state.options.mass / stepper.momentum(state.stepping));
    state.stepping.t += h * derivative;
    if (state.stepping.covTransport) {
      state.stepping.derivative(3) = derivative;
    }
  }

  template <typename propagator_state_t, typename stepper_t>
  bool transportMatrix(propagator_state_t& state, const stepper_t& stepper,
                       const double h, FreeMatrix& D) const {

    auto& sd = state.stepping.stepData;
    auto dir = stepper.direction(state.stepping);
    auto qop =
        stepper.charge(state.stepping) / stepper.momentum(state.stepping);

    D = FreeMatrix::Identity();

    double half_h = h * 0.5;
    // This sets the reference to the sub matrices
    // dFdx is already initialised as (3x3) idendity
    auto dFdT = D.block<3, 3>(0, 4);
    auto dFdL = D.block<3, 1>(0, 7);
    // dGdx is already initialised as (3x3) zero
    auto dGdT = D.block<3, 3>(4, 4);
    auto dGdL = D.block<3, 1>(4, 7);

    ActsMatrixD<3, 3> dk1dT = ActsMatrixD<3, 3>::Zero();
    ActsMatrixD<3, 3> dk2dT = ActsMatrixD<3, 3>::Identity();
    ActsMatrixD<3, 3> dk3dT = ActsMatrixD<3, 3>::Identity();
    ActsMatrixD<3, 3> dk4dT = ActsMatrixD<3, 3>::Identity();

    Vector3D dk1dL = Vector3D::Zero();
    Vector3D dk2dL = Vector3D::Zero();
    Vector3D dk3dL = Vector3D::Zero();
    Vector3D dk4dL = Vector3D::Zero();

    // For the case without energy loss
    dk1dL = dir.cross(sd.B_first);
    dk2dL = (dir + half_h * sd.k1).cross(sd.B_middle) +
            qop * half_h * dk1dL.cross(sd.B_middle);
    dk3dL = (dir + half_h * sd.k2).cross(sd.B_middle) +
            qop * half_h * dk2dL.cross(sd.B_middle);
    dk4dL =
        (dir + h * sd.k3).cross(sd.B_last) + qop * h * dk3dL.cross(sd.B_last);

    dk1dT(0, 1) = sd.B_first.z();
    dk1dT(0, 2) = -sd.B_first.y();
    dk1dT(1, 0) = -sd.B_first.z();
    dk1dT(1, 2) = sd.B_first.x();
    dk1dT(2, 0) = sd.B_first.y();
    dk1dT(2, 1) = -sd.B_first.x();
    dk1dT *= qop;

    dk2dT += half_h * dk1dT;
    dk2dT = qop * VectorHelpers::cross(dk2dT, sd.B_middle);

    dk3dT += half_h * dk2dT;
    dk3dT = qop * VectorHelpers::cross(dk3dT, sd.B_middle);

    dk4dT += h * dk3dT;
    dk4dT = qop * VectorHelpers::cross(dk4dT, sd.B_last);

    dFdT.setIdentity();
    dFdT += h / 6. * (dk1dT + dk2dT + dk3dT);
    dFdT *= h;

    dFdL = (h * h) / 6. * (dk1dL + dk2dL + dk3dL);

    dGdT += h / 6. * (dk1dT + 2. * (dk2dT + dk3dT) + dk4dT);

    dGdL = h / 6. * (dk1dL + 2. * (dk2dL + dk3dL) + dk4dL);

    D(3, 7) =
        h * state.options.mass * state.options.mass *
        stepper.charge(state.stepping) /
        (stepper.momentum(state.stepping) *
         std::hypot(1., state.options.mass / stepper.momentum(state.stepping)));
    return true;
  }
};
}  // namespace Acts