Type unrolling jacobian transformation and covariance transport

core/CMakeLists.txt

@@ -71,13 +71,13 @@ detray_add_library( detray_core core
    "include/detray/materials/material_slab.hpp"
    "include/detray/materials/predefined_materials.hpp"
    # propagator include(s)
-   "include/detray/propagator/aborters.hpp"
+   "include/detray/propagator/actors/aborters.hpp"
+   "include/detray/propagator/actors/parameter_transporter.hpp"
+   "include/detray/propagator/actors/parameter_resetter.hpp"
    "include/detray/propagator/actor_chain.hpp"
    "include/detray/propagator/base_actor.hpp"
    "include/detray/propagator/base_stepper.hpp"
    "include/detray/propagator/constrained_step.hpp"
-   "include/detray/propagator/detail/covariance_engine.hpp"
-   "include/detray/propagator/detail/jacobian_engine.hpp"
    "include/detray/propagator/line_stepper.hpp"
    "include/detray/propagator/navigation_policies.hpp"
    "include/detray/propagator/navigator.hpp"
@@ -107,7 +107,6 @@ detray_add_library( detray_core core
    "include/detray/tools/local_object_finder.hpp"
    # tracks include(s)
    "include/detray/tracks/detail/track_helper.hpp"
-   "include/detray/tracks/detail/trigonometrics.hpp"
    "include/detray/tracks/bound_track_parameters.hpp"
    "include/detray/tracks/free_track_parameters.hpp"
    "include/detray/tracks/tracks.hpp"

core/include/detray/coordinates/cartesian2.hpp

@@ -21,29 +21,31 @@ struct cartesian2 final : public coordinate_base<cartesian2, transform3_t> {
     /// @name Type definitions for the struct
     /// @{
 
-    /// Base type
+    // Base type
     using base_type = coordinate_base<cartesian2, transform3_t>;
-    /// Sclar type
+    // Sclar type
     using scalar_type = typename base_type::scalar_type;
-    /// Point in 2D space
+    // Point in 2D space
     using point2 = typename base_type::point2;
-    /// Point in 3D space
+    // Point in 3D space
     using point3 = typename base_type::point3;
-    /// Vector in 3D space
+    // Vector in 3D space
     using vector3 = typename base_type::vector3;
-    /// Matrix actor
+    // Matrix operator
     using matrix_operator = typename base_type::matrix_operator;
-    /// Matrix size type
+    // Rotation Matrix
+    using rotation_matrix = typename base_type::rotation_matrix;
+    // Matrix size type
     using size_type = typename base_type::size_type;
-    /// 2D matrix type
+    // 2D matrix type
     template <size_type ROWS, size_type COLS>
     using matrix_type = typename base_type::template matrix_type<ROWS, COLS>;
-
-    // Trigonometrics
-    using trigonometrics = typename base_type::trigonometrics;
     // Vector types
     using bound_vector = typename base_type::bound_vector;
     using free_vector = typename base_type::free_vector;
+    // Matrix types
+    using free_to_bound_matrix = typename base_type::free_to_bound_matrix;
+    using bound_to_free_matrix = typename base_type::bound_to_free_matrix;
 
     /// @}
 
@@ -63,40 +65,82 @@ struct cartesian2 final : public coordinate_base<cartesian2, transform3_t> {
     /** This method transform from a point from global cartesian 3D frame to a
      * local 2D cartesian point */
     DETRAY_HOST_DEVICE
-    inline point2 global_to_local(const transform3_t &trf, const point3 &p,
+    inline point2 global_to_local(const transform3_t &trf3, const point3 &p,
                                   const vector3 & /*d*/) const {
-        const auto local3 = trf.point_to_local(p);
+        const auto local3 = trf3.point_to_local(p);
         return this->operator()(local3);
     }
 
     /** This method transform from a local 2D cartesian point to a point global
      * cartesian 3D frame*/
     template <typename mask_t>
     DETRAY_HOST_DEVICE inline point3 local_to_global(
-        const transform3_t &trf, const mask_t & /*mask*/, const point2 &p,
+        const transform3_t &trf3, const mask_t & /*mask*/, const point2 &p,
         const vector3 & /*d*/) const {
-        return trf.point_to_global(point3{p[0], p[1], 0.});
+        return trf3.point_to_global(point3{p[0], p[1], 0.});
     }
 
-    DETRAY_HOST_DEVICE
-    inline matrix_type<3, 2> bound_to_free_rotation(
-        const transform3_t &trf3, const trigonometrics & /*t*/) {
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline vector3 normal(const transform3_t &trf3,
+                                             const mask_t & /*mask*/,
+                                             const point3 & /*pos*/,
+                                             const vector3 & /*dir*/) const {
+        vector3 ret;
+        const matrix_type<3, 1> n =
+            matrix_operator().template block<3, 1>(trf3.matrix(), 0, 2);
+        ret[0] = matrix_operator().element(n, 0, 0);
+        ret[1] = matrix_operator().element(n, 1, 0);
+        ret[2] = matrix_operator().element(n, 2, 0);
+        return ret;
+    }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline rotation_matrix reference_frame(
+        const transform3_t &trf3, const mask_t & /*mask*/,
+        const point3 & /*pos*/, const vector3 & /*dir*/) const {
+        return trf3.rotation();
+    }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_bound_pos_to_free_pos_derivative(
+        bound_to_free_matrix &bound_to_free_jacobian, const transform3_t &trf3,
+        const mask_t &mask, const point3 &pos, const vector3 &dir) const {
+
+        const rotation_matrix frame = reference_frame(trf3, mask, pos, dir);
 
         // Get d(x,y,z)/d(loc0, loc1)
-        return matrix_operator().template block<3, 2>(trf3.matrix(), 0, 0);
+        const matrix_type<3, 2> bound_pos_to_free_pos_derivative =
+            matrix_operator().template block<3, 2>(frame, 0, 0);
+
+        matrix_operator().template set_block(bound_to_free_jacobian,
+                                             bound_pos_to_free_pos_derivative,
+                                             e_free_pos0, e_bound_loc0);
     }
 
-    DETRAY_HOST_DEVICE
-    inline matrix_type<2, 3> free_to_bound_rotation(
-        const transform3_t &trf3, const trigonometrics & /*t*/) {
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_free_pos_to_bound_pos_derivative(
+        free_to_bound_matrix &free_to_bound_jacobian, const transform3_t &trf3,
+        const mask_t &mask, const point3 &pos, const vector3 &dir) const {
 
-        // Get transpose of transform3 matrix
-        const auto trf3T = matrix_operator().transpose(trf3);
+        const rotation_matrix frame = reference_frame(trf3, mask, pos, dir);
+        const rotation_matrix frameT = matrix_operator().transpose(frame);
 
         // Get d(loc0, loc1)/d(x,y,z)
-        return matrix_operator().template block<2, 3>(trf3T.matrix(), 0, 0);
+        const matrix_type<2, 3> free_pos_to_bound_pos_derivative =
+            matrix_operator().template block<2, 3>(frameT, 0, 0);
+
+        matrix_operator().template set_block(free_to_bound_jacobian,
+                                             free_pos_to_bound_pos_derivative,
+                                             e_bound_loc0, e_free_pos0);
     }
 
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_bound_angle_to_free_pos_derivative(
+        bound_to_free_matrix & /*bound_to_free_jacobian*/,
+        const transform3_t & /*trf3*/, const mask_t & /*mask*/,
+        const point3 & /*pos*/, const vector3 & /*dir*/) const {
+        // Do nothing
+    }
 };  // struct cartesian2
 
 }  // namespace detray