From 15ae220d6b699ac1057c203010f365b3e84c39f8 Mon Sep 17 00:00:00 2001 From: beomki-yeo Date: Fri, 16 Feb 2024 02:25:22 +0100 Subject: [PATCH 01/13] Backup --- .../propagator/actors}/random_scatterer.hpp | 4 +-- .../detray/propagator/base_stepper.hpp | 36 +++++++++++++++---- core/include/detray/propagator/rk_stepper.hpp | 14 +++++++- core/include/detray/propagator/rk_stepper.ipp | 35 +++++++++++++++--- .../detray/propagator}/scattering_helper.hpp | 0 .../detray/propagator/stepping_config.hpp | 3 ++ .../detray/utils}/landau_distribution.hpp | 0 tests/unit_tests/cpu/material/energy_loss.cpp | 4 +-- .../cpu/simulation/landau_sampling.cpp | 2 +- .../unit_tests/cpu/simulation/scattering.cpp | 4 +-- 10 files changed, 83 insertions(+), 19 deletions(-) rename {utils/include/detray/simulation => core/include/detray/propagator/actors}/random_scatterer.hpp (98%) rename {utils/include/detray/simulation => core/include/detray/propagator}/scattering_helper.hpp (100%) rename {utils/include/detray/simulation => core/include/detray/utils}/landau_distribution.hpp (100%) diff --git a/utils/include/detray/simulation/random_scatterer.hpp b/core/include/detray/propagator/actors/random_scatterer.hpp similarity index 98% rename from utils/include/detray/simulation/random_scatterer.hpp rename to core/include/detray/propagator/actors/random_scatterer.hpp index 9f2ef0277..c168c6c05 100644 --- a/utils/include/detray/simulation/random_scatterer.hpp +++ b/core/include/detray/propagator/actors/random_scatterer.hpp @@ -14,10 +14,10 @@ #include "detray/geometry/surface.hpp" #include "detray/materials/interaction.hpp" #include "detray/propagator/base_actor.hpp" -#include "detray/simulation/landau_distribution.hpp" -#include "detray/simulation/scattering_helper.hpp" +#include "detray/propagator/scattering_helper.hpp" #include "detray/tracks/bound_track_parameters.hpp" #include "detray/utils/axis_rotation.hpp" +#include "detray/utils/landau_distribution.hpp" #include "detray/utils/ranges.hpp" #include "detray/utils/unit_vectors.hpp" diff --git a/core/include/detray/propagator/base_stepper.hpp b/core/include/detray/propagator/base_stepper.hpp index 95596e1ed..fb56c5be5 100644 --- a/core/include/detray/propagator/base_stepper.hpp +++ b/core/include/detray/propagator/base_stepper.hpp @@ -16,6 +16,9 @@ #include "detray/propagator/stepping_config.hpp" #include "detray/tracks/tracks.hpp" +// System include(s). +#include + namespace detray { namespace stepping { @@ -140,6 +143,15 @@ class base_stepper { /// The particle pdg int _pdg = 13; // default muon +#if defined(__NO_DEVICE__) + /// Random generator + std::random_device _rd{}; + std::mt19937_64 _generator{_rd()}; + void set_seed(const uint_fast64_t sd) { + _generator.seed(sd); + } +#endif + /// Set new step constraint template DETRAY_HOST_DEVICE inline void set_constraint(scalar step_size) { @@ -153,7 +165,9 @@ class base_stepper { /// @returns access to this states step constraints DETRAY_HOST_DEVICE - inline const constraint_t &constraints() const { return _constraint; } + inline const constraint_t &constraints() const { + return _constraint; + } /// @returns access to this states step constraints DETRAY_HOST_DEVICE @@ -163,7 +177,9 @@ class base_stepper { /// @returns the navigation direction DETRAY_HOST_DEVICE - inline step::direction direction() const { return _direction; } + inline step::direction direction() const { + return _direction; + } /// Remove [all] constraints template @@ -173,19 +189,27 @@ class base_stepper { /// Set next step size DETRAY_HOST_DEVICE - inline void set_step_size(const scalar step) { _step_size = step; } + inline void set_step_size(const scalar step) { + _step_size = step; + } /// @returns the current step size of this state. DETRAY_HOST_DEVICE - inline scalar step_size() const { return _step_size; } + inline scalar step_size() const { + return _step_size; + } /// @returns this states remaining path length. DETRAY_HOST_DEVICE - inline scalar path_length() const { return _path_length; } + inline scalar path_length() const { + return _path_length; + } /// @returns the stepping inspector DETRAY_HOST - inline constexpr auto &inspector() { return _inspector; } + inline constexpr auto &inspector() { + return _inspector; + } /// Call the stepping inspector DETRAY_HOST_DEVICE diff --git a/core/include/detray/propagator/rk_stepper.hpp b/core/include/detray/propagator/rk_stepper.hpp index 68ca04bc6..c8cc4c0ec 100644 --- a/core/include/detray/propagator/rk_stepper.hpp +++ b/core/include/detray/propagator/rk_stepper.hpp @@ -72,6 +72,17 @@ class rk_stepper final const magnetic_field_t& mag_field, const detector_t& det) : base_type::state(bound_params, det), _magnetic_field(mag_field) {} +/* +#if defined(__NO_DEVICE__) + /// Random generator + std::random_device _rd{}; + std::mt19937_64 _generator{_rd()}; + + void set_seed(const uint_fast64_t sd) { + _generator.seed(sd); + } +#endif +*/ /// stepping data required for RKN4 struct { vector3 b_first{0.f, 0.f, 0.f}; @@ -95,7 +106,8 @@ class rk_stepper final /// Update the track state by Runge-Kutta-Nystrom integration. DETRAY_HOST_DEVICE - inline void advance_track(); + inline void advance_track( + const stepping::config& cfg = {}); /// Update the jacobian transport from free propagation DETRAY_HOST_DEVICE diff --git a/core/include/detray/propagator/rk_stepper.ipp b/core/include/detray/propagator/rk_stepper.ipp index ae279a4bc..b742f040f 100644 --- a/core/include/detray/propagator/rk_stepper.ipp +++ b/core/include/detray/propagator/rk_stepper.ipp @@ -6,14 +6,18 @@ */ // Project include(s). +#include "detray/definitions/detail/macros.hpp" #include "detray/geometry/detector_volume.hpp" +#include "detray/materials/interaction.hpp" +#include "detray/propagator/actors/random_scatterer.hpp" template class array_t> -DETRAY_HOST_DEVICE void -detray::rk_stepper::state::advance_track() { +DETRAY_HOST_DEVICE void detray::rk_stepper< + magnetic_field_t, transform3_t, constraint_t, policy_t, inspector_t, + array_t>::state::advance_track(const detray::stepping::config& + cfg) { const auto& sd = this->_step_data; const scalar_type h{this->_step_size}; @@ -34,7 +38,22 @@ detray::rk_stepper_mat == vacuum())) { + const auto mat = this->_mat; + if (!(mat == vacuum())) { + +#if defined(__NO_DEVICE__) + // only for simulation + if (cfg.do_scattering) { + + const scalar_type projected_scattering_angle = + interaction().compute_multiple_scattering_theta0( + h / mat.X0(), this->_pdg, this->_mass, qop, track.charge()); + + dir = random_scatterer().scatter( + dir, projected_scattering_angle, this->_generator); + } +#endif + // Reference: Eq (82) of https://doi.org/10.1016/0029-554X(81)90063-X qop = qop + h_6 * (sd.dqopds[0u] + 2.f * (sd.dqopds[1u] + sd.dqopds[2u]) + @@ -385,6 +404,12 @@ detray::rk_stepper_jac_transport = D * this->_jac_transport; } @@ -728,7 +753,7 @@ DETRAY_HOST_DEVICE bool detray::rk_stepper< } // Advance track state - stepping.advance_track(); + stepping.advance_track(cfg); // Advance jacobian transport stepping.advance_jacobian(cfg); diff --git a/utils/include/detray/simulation/scattering_helper.hpp b/core/include/detray/propagator/scattering_helper.hpp similarity index 100% rename from utils/include/detray/simulation/scattering_helper.hpp rename to core/include/detray/propagator/scattering_helper.hpp diff --git a/core/include/detray/propagator/stepping_config.hpp b/core/include/detray/propagator/stepping_config.hpp index 579f1e291..48be594b1 100644 --- a/core/include/detray/propagator/stepping_config.hpp +++ b/core/include/detray/propagator/stepping_config.hpp @@ -12,6 +12,7 @@ // System include(s). #include +#include namespace detray::stepping { @@ -34,6 +35,8 @@ struct config { scalar_t path_limit{5.f * unit::m}; /// Maximum number of Runge-Kutta step trials std::size_t max_rk_updates{10000u}; + /// Include multiple scattering -> only for simulation + bool do_scattering{false}; /// Use mean energy loss (Bethe) /// if false, most probable energy loss (Landau) will be used bool use_mean_loss{true}; diff --git a/utils/include/detray/simulation/landau_distribution.hpp b/core/include/detray/utils/landau_distribution.hpp similarity index 100% rename from utils/include/detray/simulation/landau_distribution.hpp rename to core/include/detray/utils/landau_distribution.hpp diff --git a/tests/unit_tests/cpu/material/energy_loss.cpp b/tests/unit_tests/cpu/material/energy_loss.cpp index e6c1850de..744f5d0cf 100644 --- a/tests/unit_tests/cpu/material/energy_loss.cpp +++ b/tests/unit_tests/cpu/material/energy_loss.cpp @@ -10,8 +10,8 @@ #include "detray/materials/material.hpp" #include "detray/materials/material_slab.hpp" #include "detray/materials/predefined_materials.hpp" -#include "detray/simulation/landau_distribution.hpp" -#include "detray/simulation/random_scatterer.hpp" +#include "detray/utils/landau_distribution.hpp" +#include "detray/propagator/actors/random_scatterer.hpp" #include "detray/test/types.hpp" #include "detray/utils/statistics.hpp" diff --git a/tests/unit_tests/cpu/simulation/landau_sampling.cpp b/tests/unit_tests/cpu/simulation/landau_sampling.cpp index 5734c5f8a..37d1da866 100644 --- a/tests/unit_tests/cpu/simulation/landau_sampling.cpp +++ b/tests/unit_tests/cpu/simulation/landau_sampling.cpp @@ -6,7 +6,7 @@ */ // Project include(s). -#include "detray/simulation/landau_distribution.hpp" +#include "detray/utils/landau_distribution.hpp" // GTest include(s) #include diff --git a/tests/unit_tests/cpu/simulation/scattering.cpp b/tests/unit_tests/cpu/simulation/scattering.cpp index 72040b1e5..e2f3a29d8 100644 --- a/tests/unit_tests/cpu/simulation/scattering.cpp +++ b/tests/unit_tests/cpu/simulation/scattering.cpp @@ -7,8 +7,8 @@ // Project include(s). #include "detray/propagator/actors/pointwise_material_interactor.hpp" -#include "detray/simulation/random_scatterer.hpp" -#include "detray/simulation/scattering_helper.hpp" +#include "detray/propagator/actors/random_scatterer.hpp" +#include "detray/propagator/scattering_helper.hpp" #include "detray/test/types.hpp" #include "detray/tracks/tracks.hpp" #include "detray/utils/statistics.hpp" From 6aa7d80a03c8c435de86fec2c5292c0a2ff786df Mon Sep 17 00:00:00 2001 From: beomki-yeo Date: Fri, 16 Feb 2024 16:36:47 +0100 Subject: [PATCH 02/13] Backup --- .../propagator/actors/random_scatterer.hpp | 2 +- .../detray/propagator/base_stepper.hpp | 45 +++----- .../propagator/detail/random_device.hpp | 26 +++++ .../{ => detail}/scattering_helper.hpp | 0 .../detray/propagator/line_stepper.hpp | 6 +- core/include/detray/propagator/rk_stepper.hpp | 16 +-- core/include/detray/propagator/rk_stepper.ipp | 100 ++++++++++-------- .../include/detray/test/helix_navigation.hpp | 7 +- .../detray/test/straight_line_navigation.hpp | 3 +- .../cpu/material/material_interaction.cpp | 2 +- tests/unit_tests/cpu/material/energy_loss.cpp | 2 +- .../unit_tests/cpu/simulation/scattering.cpp | 2 +- .../validation/root/rk_tolerance_comparison.C | 2 +- 13 files changed, 113 insertions(+), 100 deletions(-) create mode 100644 core/include/detray/propagator/detail/random_device.hpp rename core/include/detray/propagator/{ => detail}/scattering_helper.hpp (100%) diff --git a/core/include/detray/propagator/actors/random_scatterer.hpp b/core/include/detray/propagator/actors/random_scatterer.hpp index c168c6c05..17f78f26b 100644 --- a/core/include/detray/propagator/actors/random_scatterer.hpp +++ b/core/include/detray/propagator/actors/random_scatterer.hpp @@ -14,7 +14,7 @@ #include "detray/geometry/surface.hpp" #include "detray/materials/interaction.hpp" #include "detray/propagator/base_actor.hpp" -#include "detray/propagator/scattering_helper.hpp" +#include "detray/propagator/detail/scattering_helper.hpp" #include "detray/tracks/bound_track_parameters.hpp" #include "detray/utils/axis_rotation.hpp" #include "detray/utils/landau_distribution.hpp" diff --git a/core/include/detray/propagator/base_stepper.hpp b/core/include/detray/propagator/base_stepper.hpp index fb56c5be5..2fc6d1c63 100644 --- a/core/include/detray/propagator/base_stepper.hpp +++ b/core/include/detray/propagator/base_stepper.hpp @@ -13,16 +13,19 @@ #include "detray/geometry/surface.hpp" #include "detray/propagator/actors/parameter_resetter.hpp" #include "detray/propagator/constrained_step.hpp" +#include "detray/propagator/detail/random_device.hpp" #include "detray/propagator/stepping_config.hpp" #include "detray/tracks/tracks.hpp" -// System include(s). -#include - namespace detray { namespace stepping { +/// A void random device that does nothing. +struct void_random_device { + void set_seed(const uint_fast64_t /*sd*/) { return; } +}; + /// A void inpector that does nothing. /// /// Inspectors can be plugged in to understand the current stepper state. @@ -36,6 +39,7 @@ struct void_inspector { /// Base stepper implementation template class base_stepper { @@ -140,17 +144,12 @@ class base_stepper { /// The particle mass scalar_type _mass{105.7f * unit::MeV}; + /// The particle pdg int _pdg = 13; // default muon -#if defined(__NO_DEVICE__) - /// Random generator - std::random_device _rd{}; - std::mt19937_64 _generator{_rd()}; - void set_seed(const uint_fast64_t sd) { - _generator.seed(sd); - } -#endif + /// Random device + random_device_t rand_device; /// Set new step constraint template @@ -165,9 +164,7 @@ class base_stepper { /// @returns access to this states step constraints DETRAY_HOST_DEVICE - inline const constraint_t &constraints() const { - return _constraint; - } + inline const constraint_t &constraints() const { return _constraint; } /// @returns access to this states step constraints DETRAY_HOST_DEVICE @@ -177,9 +174,7 @@ class base_stepper { /// @returns the navigation direction DETRAY_HOST_DEVICE - inline step::direction direction() const { - return _direction; - } + inline step::direction direction() const { return _direction; } /// Remove [all] constraints template @@ -189,27 +184,19 @@ class base_stepper { /// Set next step size DETRAY_HOST_DEVICE - inline void set_step_size(const scalar step) { - _step_size = step; - } + inline void set_step_size(const scalar step) { _step_size = step; } /// @returns the current step size of this state. DETRAY_HOST_DEVICE - inline scalar step_size() const { - return _step_size; - } + inline scalar step_size() const { return _step_size; } /// @returns this states remaining path length. DETRAY_HOST_DEVICE - inline scalar path_length() const { - return _path_length; - } + inline scalar path_length() const { return _path_length; } /// @returns the stepping inspector DETRAY_HOST - inline constexpr auto &inspector() { - return _inspector; - } + inline constexpr auto &inspector() { return _inspector; } /// Call the stepping inspector DETRAY_HOST_DEVICE diff --git a/core/include/detray/propagator/detail/random_device.hpp b/core/include/detray/propagator/detail/random_device.hpp new file mode 100644 index 000000000..0fbba7c63 --- /dev/null +++ b/core/include/detray/propagator/detail/random_device.hpp @@ -0,0 +1,26 @@ +/** Detray plugins library, part of the ACTS project + * + * (c) 2024 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ + +#pragma once + +// System include(s). +#include + +namespace detray { + +namespace stepping { + +struct default_random_device { + /// Random generator + std::random_device rd{}; + std::mt19937_64 generator{rd()}; + void set_seed(const uint_fast64_t sd) { generator.seed(sd); } +}; + +} // namespace stepping + +} // namespace detray \ No newline at end of file diff --git a/core/include/detray/propagator/scattering_helper.hpp b/core/include/detray/propagator/detail/scattering_helper.hpp similarity index 100% rename from core/include/detray/propagator/scattering_helper.hpp rename to core/include/detray/propagator/detail/scattering_helper.hpp diff --git a/core/include/detray/propagator/line_stepper.hpp b/core/include/detray/propagator/line_stepper.hpp index c4b3840a5..14071c40f 100644 --- a/core/include/detray/propagator/line_stepper.hpp +++ b/core/include/detray/propagator/line_stepper.hpp @@ -12,19 +12,21 @@ #include "detray/definitions/detail/qualifiers.hpp" #include "detray/navigation/policies.hpp" #include "detray/propagator/base_stepper.hpp" +#include "detray/propagator/detail/random_device.hpp" namespace detray { /// Straight line stepper implementation template class line_stepper final : public base_stepper { public: - using base_type = - base_stepper; + using base_type = base_stepper; using free_track_parameters_type = typename base_type::free_track_parameters_type; diff --git a/core/include/detray/propagator/rk_stepper.hpp b/core/include/detray/propagator/rk_stepper.hpp index c8cc4c0ec..4f2fb76ce 100644 --- a/core/include/detray/propagator/rk_stepper.hpp +++ b/core/include/detray/propagator/rk_stepper.hpp @@ -28,14 +28,15 @@ namespace detray { template class array_t = darray> class rk_stepper final : public base_stepper { public: - using base_type = - base_stepper; + using base_type = base_stepper; using transform3_type = transform3_t; using scalar_type = typename transform3_type::scalar_type; @@ -72,17 +73,6 @@ class rk_stepper final const magnetic_field_t& mag_field, const detector_t& det) : base_type::state(bound_params, det), _magnetic_field(mag_field) {} -/* -#if defined(__NO_DEVICE__) - /// Random generator - std::random_device _rd{}; - std::mt19937_64 _generator{_rd()}; - - void set_seed(const uint_fast64_t sd) { - _generator.seed(sd); - } -#endif -*/ /// stepping data required for RKN4 struct { vector3 b_first{0.f, 0.f, 0.f}; diff --git a/core/include/detray/propagator/rk_stepper.ipp b/core/include/detray/propagator/rk_stepper.ipp index b742f040f..6522031f3 100644 --- a/core/include/detray/propagator/rk_stepper.ipp +++ b/core/include/detray/propagator/rk_stepper.ipp @@ -12,12 +12,12 @@ #include "detray/propagator/actors/random_scatterer.hpp" template class array_t> -DETRAY_HOST_DEVICE void detray::rk_stepper< - magnetic_field_t, transform3_t, constraint_t, policy_t, inspector_t, - array_t>::state::advance_track(const detray::stepping::config& - cfg) { + typename constraint_t, typename policy_t, typename random_device_t, + typename inspector_t, template class array_t> +DETRAY_HOST_DEVICE void +detray::rk_stepper::state:: + advance_track(const detray::stepping::config& cfg) { const auto& sd = this->_step_data; const scalar_type h{this->_step_size}; @@ -41,18 +41,22 @@ DETRAY_HOST_DEVICE void detray::rk_stepper< const auto mat = this->_mat; if (!(mat == vacuum())) { -#if defined(__NO_DEVICE__) - // only for simulation - if (cfg.do_scattering) { + // only for simulation with valid random device + if constexpr (!std::is_same::value) { - const scalar_type projected_scattering_angle = - interaction().compute_multiple_scattering_theta0( - h / mat.X0(), this->_pdg, this->_mass, qop, track.charge()); + if (cfg.do_scattering) { + const scalar_type projected_scattering_angle = + interaction() + .compute_multiple_scattering_theta0( + h / mat.X0(), this->_pdg, this->_mass, qop, + track.charge()); - dir = random_scatterer().scatter( - dir, projected_scattering_angle, this->_generator); + dir = random_scatterer().scatter( + dir, projected_scattering_angle, + this->rand_device.generator); + } } -#endif // Reference: Eq (82) of https://doi.org/10.1016/0029-554X(81)90063-X qop = @@ -67,11 +71,11 @@ DETRAY_HOST_DEVICE void detray::rk_stepper< } template class array_t> + typename constraint_t, typename policy_t, typename random_device_t, + typename inspector_t, template class array_t> DETRAY_HOST_DEVICE void detray::rk_stepper::state:: + random_device_t, inspector_t, array_t>::state:: advance_jacobian(const detray::stepping::config& cfg) { /// The calculations are based on ATL-SOFT-PUB-2009-002. The update of the /// Jacobian matrix is requires only the calculation of eq. 17 and 18. @@ -414,11 +418,11 @@ detray::rk_stepper class array_t> + typename constraint_t, typename policy_t, typename random_device_t, + typename inspector_t, template class array_t> DETRAY_HOST_DEVICE auto detray::rk_stepper::state:: + random_device_t, inspector_t, array_t>::state:: evaluate_dqopds(const std::size_t i, const typename transform3_t::scalar_type h, const scalar dqopds_prev, @@ -453,10 +457,11 @@ detray::rk_stepper class array_t> + typename constraint_t, typename policy_t, typename random_device_t, + typename inspector_t, template class array_t> DETRAY_HOST_DEVICE auto detray::rk_stepper< - magnetic_field_t, transform3_t, constraint_t, policy_t, inspector_t, + magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t, + inspector_t, array_t>::state::evaluate_dtds(const vector3& b_field, const std::size_t i, const typename transform3_t::scalar_type h, const vector3& dtds_prev, @@ -478,11 +483,11 @@ DETRAY_HOST_DEVICE auto detray::rk_stepper< } template class array_t> + typename constraint_t, typename policy_t, typename random_device_t, + typename inspector_t, template class array_t> DETRAY_HOST_DEVICE auto detray::rk_stepper< - magnetic_field_t, transform3_t, constraint_t, policy_t, inspector_t, - array_t>::state::evaluate_field_gradient(const vector3& pos) + magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t, + inspector_t, array_t>::state::evaluate_field_gradient(const vector3& pos) -> matrix_type<3, 3> { matrix_type<3, 3> dBdr = matrix_operator().template zero<3, 3>(); @@ -520,22 +525,22 @@ DETRAY_HOST_DEVICE auto detray::rk_stepper< } template class array_t> + typename constraint_t, typename policy_t, typename random_device_t, + typename inspector_t, template class array_t> DETRAY_HOST_DEVICE auto detray::rk_stepper::state::dqopds() const -> - typename transform3_t::scalar_type { + random_device_t, inspector_t, array_t>::state::dqopds() const + -> typename transform3_t::scalar_type { return this->_step_data.dqopds[3u]; } template class array_t> -DETRAY_HOST_DEVICE auto -detray::rk_stepper::state::dqopds(const scalar_type qop) - const -> typename transform3_t::scalar_type { + typename constraint_t, typename policy_t, typename random_device_t, + typename inspector_t, template class array_t> +DETRAY_HOST_DEVICE auto detray::rk_stepper< + magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t, + inspector_t, array_t>::state::dqopds(const scalar_type qop) const -> + typename transform3_t::scalar_type { const auto& mat = this->_mat; @@ -564,11 +569,11 @@ detray::rk_stepper class array_t> + typename constraint_t, typename policy_t, typename random_device_t, + typename inspector_t, template class array_t> DETRAY_HOST_DEVICE auto detray::rk_stepper< - magnetic_field_t, transform3_t, constraint_t, policy_t, inspector_t, - array_t>::state::d2qopdsdqop(const scalar_type qop) const -> + magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t, + inspector_t, array_t>::state::d2qopdsdqop(const scalar_type qop) const -> typename transform3_t::scalar_type { using scalar_t = typename transform3_t::scalar_type; @@ -607,13 +612,14 @@ DETRAY_HOST_DEVICE auto detray::rk_stepper< } template class array_t> + typename constraint_t, typename policy_t, typename random_device_t, + typename inspector_t, template class array_t> template DETRAY_HOST_DEVICE bool detray::rk_stepper< - magnetic_field_t, transform3_t, constraint_t, policy_t, inspector_t, - array_t>::step(propagation_state_t& propagation, - const detray::stepping::config& cfg) { + magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t, + inspector_t, array_t>::step(propagation_state_t& propagation, + const detray::stepping::config& + cfg) { // Get stepper and navigator states state& stepping = propagation._stepping; diff --git a/tests/include/detray/test/helix_navigation.hpp b/tests/include/detray/test/helix_navigation.hpp index bf8c1d7f8..3c785c55b 100644 --- a/tests/include/detray/test/helix_navigation.hpp +++ b/tests/include/detray/test/helix_navigation.hpp @@ -107,9 +107,10 @@ class helix_navigation : public test::fixture_base<> { using navigator_t = navigator; // Runge-Kutta stepper using bfield_t = bfield::const_field_t; - using stepper_t = rk_stepper; + using stepper_t = + rk_stepper; // Propagator with pathlimit aborter using actor_chain_t = actor_chain; using propagator_t = propagator; diff --git a/tests/include/detray/test/straight_line_navigation.hpp b/tests/include/detray/test/straight_line_navigation.hpp index dc59d200b..d49c77fc8 100644 --- a/tests/include/detray/test/straight_line_navigation.hpp +++ b/tests/include/detray/test/straight_line_navigation.hpp @@ -110,7 +110,8 @@ class straight_line_navigation : public test::fixture_base<> { // Line stepper using stepper_t = line_stepper; + stepper_default_policy, stepping::void_random_device, + stepping::print_inspector>; // Propagator with pathlimit aborter using actor_chain_t = actor_chain; using propagator_t = propagator; diff --git a/tests/integration_tests/cpu/material/material_interaction.cpp b/tests/integration_tests/cpu/material/material_interaction.cpp index f5ba6f374..b2252cd36 100644 --- a/tests/integration_tests/cpu/material/material_interaction.cpp +++ b/tests/integration_tests/cpu/material/material_interaction.cpp @@ -22,10 +22,10 @@ #include "detray/propagator/actors/parameter_resetter.hpp" #include "detray/propagator/actors/parameter_transporter.hpp" #include "detray/propagator/actors/pointwise_material_interactor.hpp" +#include "detray/propagator/actors/random_scatterer.hpp" #include "detray/propagator/line_stepper.hpp" #include "detray/propagator/propagator.hpp" #include "detray/propagator/rk_stepper.hpp" -#include "detray/simulation/random_scatterer.hpp" #include "detray/test/types.hpp" #include "detray/utils/inspectors.hpp" #include "detray/utils/statistics.hpp" diff --git a/tests/unit_tests/cpu/material/energy_loss.cpp b/tests/unit_tests/cpu/material/energy_loss.cpp index 744f5d0cf..9b644298e 100644 --- a/tests/unit_tests/cpu/material/energy_loss.cpp +++ b/tests/unit_tests/cpu/material/energy_loss.cpp @@ -10,9 +10,9 @@ #include "detray/materials/material.hpp" #include "detray/materials/material_slab.hpp" #include "detray/materials/predefined_materials.hpp" -#include "detray/utils/landau_distribution.hpp" #include "detray/propagator/actors/random_scatterer.hpp" #include "detray/test/types.hpp" +#include "detray/utils/landau_distribution.hpp" #include "detray/utils/statistics.hpp" // GTest include(s). diff --git a/tests/unit_tests/cpu/simulation/scattering.cpp b/tests/unit_tests/cpu/simulation/scattering.cpp index e2f3a29d8..912aa7928 100644 --- a/tests/unit_tests/cpu/simulation/scattering.cpp +++ b/tests/unit_tests/cpu/simulation/scattering.cpp @@ -8,7 +8,7 @@ // Project include(s). #include "detray/propagator/actors/pointwise_material_interactor.hpp" #include "detray/propagator/actors/random_scatterer.hpp" -#include "detray/propagator/scattering_helper.hpp" +#include "detray/propagator/detail/scattering_helper.hpp" #include "detray/test/types.hpp" #include "detray/tracks/tracks.hpp" #include "detray/utils/statistics.hpp" diff --git a/tests/validation/root/rk_tolerance_comparison.C b/tests/validation/root/rk_tolerance_comparison.C index b6709de03..d5744f4d4 100644 --- a/tests/validation/root/rk_tolerance_comparison.C +++ b/tests/validation/root/rk_tolerance_comparison.C @@ -217,7 +217,7 @@ void draw_mean_step_size(const std::string header_title, gr->GetXaxis()->SetTitle("log_{10}(#tau [mm])"); gr->GetYaxis()->SetTitle("log_{10}(Mean of average step sizes [mm])"); gr->GetXaxis()->SetLimits(x_vec.front() - 0.5, x_vec.back() + 0.5); - gr->GetYaxis()->SetRangeUser(means.front() - 0.2, means.back() + 0.95); + gr->GetYaxis()->SetRangeUser(means.front() - 0.2, means.back() + 1.00); gr->GetXaxis()->SetLabelSize(label_font_size); gr->GetYaxis()->SetLabelSize(label_font_size); gr->GetXaxis()->SetTitleSize(title_font_size); From 5c6f2dea22d7f19109ff9b06cec86a9f14a25a17 Mon Sep 17 00:00:00 2001 From: beomki-yeo Date: Fri, 16 Feb 2024 19:05:34 +0100 Subject: [PATCH 03/13] Backup --- tests/scripts/run_jacobian_validation.sh | 16 ++++++++++++---- tests/validation/root/covariance_validation.C | 11 ++++++----- 2 files changed, 18 insertions(+), 9 deletions(-) diff --git a/tests/scripts/run_jacobian_validation.sh b/tests/scripts/run_jacobian_validation.sh index d3542471d..ed362cef8 100644 --- a/tests/scripts/run_jacobian_validation.sh +++ b/tests/scripts/run_jacobian_validation.sh @@ -16,6 +16,9 @@ n_tracks_per_thread=1000 log10_min_rk_tol=-6 log10_max_rk_tol=2 +# log10 rk tolerance for covariance transport +log10_rk_tol=-4 + # Helix intersector tolerance [log10 in mm] log10_helix_tol=-3 # Surface tolerance [log10 in mm] @@ -30,7 +33,7 @@ skip_second_phase=false # Verbose level verbose_level=1 -while getopts "hd:n:t:p:q:i:s:r:v:" arg; do +while getopts "hd:n:t:p:q:c:i:s:r:v:" arg; do case $arg in h) echo "" @@ -42,6 +45,7 @@ while getopts "hd:n:t:p:q:i:s:r:v:" arg; do echo "-t " echo "-p " echo "-q " + echo "-c " echo "-i " echo "-s " echo "-r " @@ -69,6 +73,10 @@ while getopts "hd:n:t:p:q:i:s:r:v:" arg; do log10_max_rk_tol=$OPTARG echo "log10(max_rk_error_tolerance_in_mm): ${log10_max_rk_tol}" ;; + c) + log10_rk_tol=$OPTARG + echo "log10(rk_error_tolerance_in_mm_for_covariance_transport): ${log10_rk_tol}" + ;; i) log10_helix_tol=$OPTARG log10_on_surface_tol=$OPTARG @@ -147,7 +155,7 @@ if [ "$skip_second_phase" = false ] ; then --rk-tolerance-iterate-mode=false \ --n-tracks=${n_tracks_per_thread} \ --n-skips=${n_skips} \ - --log10-rk-tolerance=${log10_min_rk_tol} \ + --log10-rk-tolerance=${log10_rk_tol} \ --log10-helix-tolerance=${log10_helix_tol} \ --log10-on-surface-tolerance=${log10_on_surface_tol} \ --mc-seed=${mc_seed} @@ -211,8 +219,8 @@ if [ "$skip_second_phase" = false ] ; then # Run covariance_validation.C root -q -l ../../../tests/validation/root/covariance_validation.C+O -else - +else + # Run rk_tolerance_comparision.C root '../../../tests/validation/root/rk_tolerance_comparison.C+O('${log10_min_rk_tol}','${log10_max_rk_tol}')' fi diff --git a/tests/validation/root/covariance_validation.C b/tests/validation/root/covariance_validation.C index 4cb858bda..5c7b5d9c3 100644 --- a/tests/validation/root/covariance_validation.C +++ b/tests/validation/root/covariance_validation.C @@ -30,6 +30,9 @@ #include namespace { +double pull_min = -8.; +double pull_max = 8.; +double bin_width = 0.1; double x_pos = 0.16f; double title_x = x_pos; double title_y = 0.8f; @@ -124,8 +127,8 @@ std::pair fit_pval(TH1D* h_pval) { } void set_yaxis_title(TH1D* h, const double text_size) { - double bin_width = h->GetBinWidth(0u); - std::string str = std::to_string(bin_width); + double bw = h->GetBinWidth(0u); + std::string str = std::to_string(bw); str.erase(str.find_last_not_of('0') + 1, std::string::npos); str.erase(str.find_last_not_of('.') + 1, std::string::npos); std::string y_axis_title = "Counts / (" + str + ")"; @@ -322,9 +325,7 @@ void read_tree(TTree* t, const std::string& tag, const std::string& title) { const std::array cdim1{700, 500}; const std::array cdim2{700, 500}; - float pull_min = -5.f; - float pull_max = 5.f; - int n_bins = 100; + int n_bins = (pull_max - pull_min) / bin_width; double pull_l0; double pull_l1; From 2f5aa683513c9c9cab472d604643f4e0035d14f2 Mon Sep 17 00:00:00 2001 From: beomki-yeo Date: Sat, 17 Feb 2024 00:21:33 +0100 Subject: [PATCH 04/13] Backup --- core/include/detray/propagator/rk_stepper.ipp | 6 -- core/include/detray/utils/axis_rotation.hpp | 2 + core/include/detray/utils/matrix_helper.hpp | 97 +++++++++++++++++++ tests/unit_tests/cpu/utils/matrix_helper.cpp | 45 ++++++++- 4 files changed, 143 insertions(+), 7 deletions(-) diff --git a/core/include/detray/propagator/rk_stepper.ipp b/core/include/detray/propagator/rk_stepper.ipp index 6522031f3..7c4d80a13 100644 --- a/core/include/detray/propagator/rk_stepper.ipp +++ b/core/include/detray/propagator/rk_stepper.ipp @@ -408,12 +408,6 @@ detray::rk_stepper_jac_transport = D * this->_jac_transport; } diff --git a/core/include/detray/utils/axis_rotation.hpp b/core/include/detray/utils/axis_rotation.hpp index 5e40817ce..c45b24373 100644 --- a/core/include/detray/utils/axis_rotation.hpp +++ b/core/include/detray/utils/axis_rotation.hpp @@ -41,6 +41,8 @@ struct axis_rotation { scalar_type cos_theta{math::cos(theta)}; + // Rodrigues' rotation_formula + // https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula matrix_type<3, 3> I = matrix_operator().template identity<3, 3>(); matrix_type<3, 3> axis_cross = mat_helper().cross_matrix(U); matrix_type<3, 3> axis_outer = mat_helper().outer_product(U, U); diff --git a/core/include/detray/utils/matrix_helper.hpp b/core/include/detray/utils/matrix_helper.hpp index 5038ab01e..b7c3c787c 100644 --- a/core/include/detray/utils/matrix_helper.hpp +++ b/core/include/detray/utils/matrix_helper.hpp @@ -67,6 +67,12 @@ struct matrix_helper { } /// Cross product matrix + /// + /// For a given vector (v) = (x,y,z), return the following matrix: + /// + /// [ 0 -z y ] + /// [ z 0 -x ] + /// [-y x 0 ] DETRAY_HOST_DEVICE inline matrix_type<3, 3> cross_matrix(const vector3& v) const { matrix_type<3, 3> ret; @@ -84,6 +90,8 @@ struct matrix_helper { } /// Outer product operation + /// + /// @return 3x3 matrix M = v * v^T where v is a 3-elements column matrix DETRAY_HOST_DEVICE inline matrix_type<3, 3> outer_product(const vector3& v1, const vector3& v2) const { @@ -101,6 +109,8 @@ struct matrix_helper { } /// Cholesky decompose + /// + /// @return the lower triangle matrix (L) that satisfies M = L * L^T template DETRAY_HOST_DEVICE inline matrix_type cholesky_decompose( const matrix_type& mat) const { @@ -127,6 +137,93 @@ struct matrix_helper { return L; } + + /// @return the rotation matrix (R) that converts the local curvilinear + /// cartesian to the global cartesian + /// + /// The local z-axis (w) is defined as the given vector. If we say that u + /// and v are the local x and y axis, R is defined as follows: + /// + /// [ u₁ v₁ w₁ ] + /// R = [ u₂ v₂ w₂ ] + /// [ u₃ v₃ w₃ ] + /// + /// One way to obtain R is solving Rodrigues' rotation_formula (This is the + /// formula implemented in axis_rotation.hpp): + /// + /// cos_α I + sin_α cross_matrix(k) + (1 - cos_α) * outer_product(k) + /// + /// k: normalization of (0,0,1) X w + /// α: the angle required to rotate (0,0,1) to w around the axis of k + /// cross_matrix: function defined above + /// outer_product: function defined_above + /// + /// + /// In a nutshell, the rest of R elements is calculated as: + /// + /// u₁ = (w₁ w₁ w₃ + w₂ w₂) / ( 1 - w₃ w₃ ) + /// u₂ = -w₁w₂/(1 + w₃) + /// u₃ = -w₁ + /// v₁ = -w₁w₂/(1 + w₃) + /// v₂ = (w₁ w₁ + w₂ w₂ w₃) / ( 1 - w₃ w₃ ) + /// v₃ = -w₂ + DETRAY_HOST_DEVICE + inline matrix_type<3, 3> curvilinear_to_global(const vector3& w) const { + + matrix_type<3, 3> R = matrix_operator().template zero<3, 3>(); + + const scalar_type w1w2 = w[0] * w[1]; + + // Set u + getter::element(R, 0u, 0u) = + (w[0] * w[0] * w[2] + w[1] * w[1]) / (1 - w[2] * w[2]); + getter::element(R, 1u, 0u) = -w1w2 / (1 + w[2]); + getter::element(R, 2u, 0u) = -w[0]; + + // Set v + getter::element(R, 0u, 1u) = getter::element(R, 1u, 0u); + getter::element(R, 1u, 1u) = + (w[0] * w[0] + w[1] * w[1] * w[2]) / (1 - w[2] * w[2]); + getter::element(R, 2u, 1u) = -w[1]; + + // Set w + getter::element(R, 0u, 2u) = w[0u]; + getter::element(R, 1u, 2u) = w[1u]; + getter::element(R, 2u, 2u) = w[2u]; + + return R; + } + + /// @return the rotation matrix (R) that converts the global cartesian to + /// the local curvilinear cartesian + /// + /// It is the transpose of curvlinear_to_global + DETRAY_HOST_DEVICE + inline matrix_type<3, 3> global_to_curvilinear(const vector3& w) { + + matrix_type<3, 3> R = matrix_operator().template zero<3, 3>(); + + const scalar_type w1w2 = w[0] * w[1]; + + // Set u + getter::element(R, 0u, 0u) = + (w[0] * w[0] * w[2] + w[1] * w[1]) / (1 - w[2] * w[2]); + getter::element(R, 1u, 0u) = -w1w2 / (1 + w[2]); + getter::element(R, 2u, 0u) = w[0]; + + // Set v + getter::element(R, 0u, 1u) = getter::element(R, 1u, 0u); + getter::element(R, 1u, 1u) = + (w[0] * w[0] + w[1] * w[1] * w[2]) / (1 - w[2] * w[2]); + getter::element(R, 2u, 1u) = w[1]; + + // Set w + getter::element(R, 0u, 2u) = -w[0u]; + getter::element(R, 1u, 2u) = -w[1u]; + getter::element(R, 2u, 2u) = w[2u]; + + return R; + } }; } // namespace detray diff --git a/tests/unit_tests/cpu/utils/matrix_helper.cpp b/tests/unit_tests/cpu/utils/matrix_helper.cpp index c749523ec..71862987c 100644 --- a/tests/unit_tests/cpu/utils/matrix_helper.cpp +++ b/tests/unit_tests/cpu/utils/matrix_helper.cpp @@ -165,4 +165,47 @@ GTEST_TEST(detray_utils, cholesky_decomposition) { static_cast(getter::element(B, i, j))); } } -} \ No newline at end of file +} + +GTEST_TEST(detray_utils, curvilinear) { + + const vector3 t1 = vector::normalize(vector3{1.f, 2.f, 3.f}); + const vector3 z{0.f, 0.f, 1.f}; + + // Convert the (0,0,1) in local coordinate to global cooridnate + const matrix_type<3, 3> R1 = + matrix_helper().curvilinear_to_global(t1); + const vector3 t2 = R1 * z; + + // The converted vector should be eqaul to the local-z axis in global + // coordinate + EXPECT_NEAR(static_cast(t1[0]), static_cast(t2[0]), + tolerance); + EXPECT_NEAR(static_cast(t1[1]), static_cast(t2[1]), + tolerance); + EXPECT_NEAR(static_cast(t1[2]), static_cast(t2[2]), + tolerance); + + // Convert the local-z axis in global coordinate to local coordinate + const matrix_type<3, 3> R2 = + matrix_helper().global_to_curvilinear(t1); + const vector3 t3 = R2 * t1; + + // The converted vector should be equal to unit z vector + EXPECT_NEAR(static_cast(z[0]), static_cast(t3[0]), tolerance); + EXPECT_NEAR(static_cast(z[1]), static_cast(t3[1]), tolerance); + EXPECT_NEAR(static_cast(z[2]), static_cast(t3[2]), tolerance); + + // Round trip test + const matrix_type<3, 3> I = R1 * R2; + + EXPECT_NEAR(static_cast(getter::element(I, 0, 0)), 1.f, tolerance); + EXPECT_NEAR(static_cast(getter::element(I, 0, 1)), 0.f, tolerance); + EXPECT_NEAR(static_cast(getter::element(I, 0, 2)), 0.f, tolerance); + EXPECT_NEAR(static_cast(getter::element(I, 1, 0)), 0.f, tolerance); + EXPECT_NEAR(static_cast(getter::element(I, 1, 1)), 1.f, tolerance); + EXPECT_NEAR(static_cast(getter::element(I, 1, 2)), 0.f, tolerance); + EXPECT_NEAR(static_cast(getter::element(I, 2, 0)), 0.f, tolerance); + EXPECT_NEAR(static_cast(getter::element(I, 2, 1)), 0.f, tolerance); + EXPECT_NEAR(static_cast(getter::element(I, 2, 2)), 1.f, tolerance); +} From 0f0c0e998cdf3f724d3094a3f0aaedb646bf49f3 Mon Sep 17 00:00:00 2001 From: beomki-yeo Date: Sat, 17 Feb 2024 00:23:47 +0100 Subject: [PATCH 05/13] Backup --- .gitignore | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.gitignore b/.gitignore index 917a2596f..e72c23650 100644 --- a/.gitignore +++ b/.gitignore @@ -10,3 +10,8 @@ # python files **/__pycache__ + +# root files +**/*.d +**/*.so +**/*C_ACLiC_dict* \ No newline at end of file From e94a71df35853b8ad1ef32367c1d617bc9365918 Mon Sep 17 00:00:00 2001 From: beomki-yeo Date: Sat, 17 Feb 2024 01:08:41 +0100 Subject: [PATCH 06/13] Backup --- .../detray/propagator/actors/parameter_transporter.hpp | 3 +++ core/include/detray/propagator/rk_stepper.hpp | 5 +++++ core/include/detray/propagator/rk_stepper.ipp | 9 +++++++++ tests/validation/root/covariance_validation.C | 4 ++-- 4 files changed, 19 insertions(+), 2 deletions(-) diff --git a/core/include/detray/propagator/actors/parameter_transporter.hpp b/core/include/detray/propagator/actors/parameter_transporter.hpp index 4f6acbf51..d19def4b9 100644 --- a/core/include/detray/propagator/actors/parameter_transporter.hpp +++ b/core/include/detray/propagator/actors/parameter_transporter.hpp @@ -119,6 +119,9 @@ struct parameter_transporter : actor { matrix_operator().transpose(stepping._full_jacobian); } + // Update covariance with the multiple scattering + stepping.add_multiple_scattering_covariance(); + // Calculate surface-to-surface covariance transport stepping._bound_params.set_covariance(new_cov); } diff --git a/core/include/detray/propagator/rk_stepper.hpp b/core/include/detray/propagator/rk_stepper.hpp index 4f2fb76ce..7817515f8 100644 --- a/core/include/detray/propagator/rk_stepper.hpp +++ b/core/include/detray/propagator/rk_stepper.hpp @@ -104,6 +104,11 @@ class rk_stepper final inline void advance_jacobian( const stepping::config& cfg = {}); + /// Update the covariance from multiple scattering + DETRAY_HOST_DEVICE + inline void add_multiple_scattering_covariance( + const stepping::config& cfg = {}); + /// evaulate dqopds for a given step size and material DETRAY_HOST_DEVICE inline scalar_type evaluate_dqopds( diff --git a/core/include/detray/propagator/rk_stepper.ipp b/core/include/detray/propagator/rk_stepper.ipp index 7c4d80a13..5451b8b0a 100644 --- a/core/include/detray/propagator/rk_stepper.ipp +++ b/core/include/detray/propagator/rk_stepper.ipp @@ -411,6 +411,15 @@ detray::rk_stepper_jac_transport = D * this->_jac_transport; } +template class array_t> +DETRAY_HOST_DEVICE void +detray::rk_stepper::state:: + add_multiple_scattering_covariance( + const detray::stepping::config& cfg) {} + template class array_t> diff --git a/tests/validation/root/covariance_validation.C b/tests/validation/root/covariance_validation.C index 5c7b5d9c3..6481d09b8 100644 --- a/tests/validation/root/covariance_validation.C +++ b/tests/validation/root/covariance_validation.C @@ -30,8 +30,8 @@ #include namespace { -double pull_min = -8.; -double pull_max = 8.; +double pull_min = -5.; +double pull_max = 5.; double bin_width = 0.1; double x_pos = 0.16f; double title_x = x_pos; From 21c938fc14c84c89baf05942e1691027325aca92 Mon Sep 17 00:00:00 2001 From: beomki-yeo Date: Sat, 17 Feb 2024 14:20:24 +0100 Subject: [PATCH 07/13] Backup --- .../detray/propagator/actors/aborters.hpp | 2 +- .../propagator/actors/parameter_resetter.hpp | 2 +- .../actors/parameter_transporter.hpp | 4 +- .../detray/propagator/base_stepper.hpp | 2 +- .../detray/propagator/line_stepper.hpp | 5 +- core/include/detray/propagator/rk_stepper.hpp | 7 +- core/include/detray/propagator/rk_stepper.ipp | 112 ++++++++++++------ .../detray/propagator/stepping_config.hpp | 6 +- .../cpu/propagator/propagator.cpp | 11 +- 9 files changed, 99 insertions(+), 52 deletions(-) diff --git a/core/include/detray/propagator/actors/aborters.hpp b/core/include/detray/propagator/actors/aborters.hpp index d478f8b70..0b07d15b6 100644 --- a/core/include/detray/propagator/actors/aborters.hpp +++ b/core/include/detray/propagator/actors/aborters.hpp @@ -109,7 +109,7 @@ struct next_surface_aborter : actor { // Abort at the next sensitive surface if (navigation.is_on_sensitive() && - stepping._s > abrt_state.min_step_length) { + stepping._path_length_per_surface > abrt_state.min_step_length) { prop_state._heartbeat &= navigation.abort(); abrt_state.success = true; } diff --git a/core/include/detray/propagator/actors/parameter_resetter.hpp b/core/include/detray/propagator/actors/parameter_resetter.hpp index 7decf874f..b9214ffc9 100644 --- a/core/include/detray/propagator/actors/parameter_resetter.hpp +++ b/core/include/detray/propagator/actors/parameter_resetter.hpp @@ -49,7 +49,7 @@ struct parameter_resetter : actor { trf3, mask, stepping._bound_params.vector())); // Reset the path length - stepping._s = 0; + stepping._path_length_per_surface = 0; // Reset jacobian coordinate transformation at the current surface stepping._jac_to_global = local_coordinate.bound_to_free_jacobian( diff --git a/core/include/detray/propagator/actors/parameter_transporter.hpp b/core/include/detray/propagator/actors/parameter_transporter.hpp index d19def4b9..91921b946 100644 --- a/core/include/detray/propagator/actors/parameter_transporter.hpp +++ b/core/include/detray/propagator/actors/parameter_transporter.hpp @@ -119,9 +119,11 @@ struct parameter_transporter : actor { matrix_operator().transpose(stepping._full_jacobian); } + /* // Update covariance with the multiple scattering stepping.add_multiple_scattering_covariance(); - + */ + // Calculate surface-to-surface covariance transport stepping._bound_params.set_covariance(new_cov); } diff --git a/core/include/detray/propagator/base_stepper.hpp b/core/include/detray/propagator/base_stepper.hpp index 2fc6d1c63..e9f9a364a 100644 --- a/core/include/detray/propagator/base_stepper.hpp +++ b/core/include/detray/propagator/base_stepper.hpp @@ -137,7 +137,7 @@ class base_stepper { /// Track path length from the last surface. It will be reset to 0 when /// the track reaches a new surface - scalar _s{0.}; + scalar _path_length_per_surface{0.}; /// Current step size scalar _step_size{0.}; diff --git a/core/include/detray/propagator/line_stepper.hpp b/core/include/detray/propagator/line_stepper.hpp index 14071c40f..1b0581cb2 100644 --- a/core/include/detray/propagator/line_stepper.hpp +++ b/core/include/detray/propagator/line_stepper.hpp @@ -58,7 +58,7 @@ class line_stepper final track.set_pos(track.pos() + track.dir() * this->_step_size); this->_path_length += this->_step_size; - this->_s += this->_step_size; + this->_path_length_per_surface += this->_step_size; } DETRAY_HOST_DEVICE @@ -80,6 +80,9 @@ class line_stepper final this->_jac_transport = D * this->_jac_transport; } + DETRAY_HOST_DEVICE + inline void add_multiple_scattering_covariance() { return; } + DETRAY_HOST_DEVICE inline vector3 dtds() const { return vector3{0.f, 0.f, 0.f}; } diff --git a/core/include/detray/propagator/rk_stepper.hpp b/core/include/detray/propagator/rk_stepper.hpp index 7817515f8..890cf20a1 100644 --- a/core/include/detray/propagator/rk_stepper.hpp +++ b/core/include/detray/propagator/rk_stepper.hpp @@ -54,6 +54,7 @@ class rk_stepper final template using matrix_type = typename matrix_operator::template matrix_type; + using interaction_type = interaction; DETRAY_HOST_DEVICE rk_stepper() {} @@ -106,15 +107,13 @@ class rk_stepper final /// Update the covariance from multiple scattering DETRAY_HOST_DEVICE - inline void add_multiple_scattering_covariance( - const stepping::config& cfg = {}); + inline void add_multiple_scattering_covariance(); /// evaulate dqopds for a given step size and material DETRAY_HOST_DEVICE inline scalar_type evaluate_dqopds( const std::size_t i, const typename transform3_t::scalar_type h, - const scalar dqopds_prev, - const detray::stepping::config& cfg); + const scalar dqopds_prev); /// evaulate dtds for runge kutta stepping DETRAY_HOST_DEVICE diff --git a/core/include/detray/propagator/rk_stepper.ipp b/core/include/detray/propagator/rk_stepper.ipp index 5451b8b0a..774c20e3c 100644 --- a/core/include/detray/propagator/rk_stepper.ipp +++ b/core/include/detray/propagator/rk_stepper.ipp @@ -45,12 +45,12 @@ detray::rk_stepper::value) { - if (cfg.do_scattering) { + // Scatter if it is for the simulation + if (cfg.is_simulation) { const scalar_type projected_scattering_angle = - interaction() - .compute_multiple_scattering_theta0( - h / mat.X0(), this->_pdg, this->_mass, qop, - track.charge()); + interaction_type().compute_multiple_scattering_theta0( + h / mat.X0(), this->_pdg, this->_mass, qop, + track.charge()); dir = random_scatterer().scatter( dir, projected_scattering_angle, @@ -59,15 +59,39 @@ detray::rk_stepper::value) { + // Apply mean probable energy loss instead + const scalar_type e_loss_mpv = + interaction_type().compute_energy_loss_landau( + this->_path_length_per_surface, mat, this->_pdg, + this->_mass, qop, track.charge()); + + const scalar_type e_loss_sigma = + interaction_type().compute_energy_loss_landau_sigma( + this->_path_length_per_surface, mat, this->_pdg, + this->_mass, qop, track.charge()); + + // Get the new momentum + const auto new_mom = attenuate(e_loss_mpv, e_loss_sigma, + this->mass, track.charge() / qop, + this->rand_device.generator); + + qop = track.charge() / new_mom; + } + } } track.set_qop(qop); // Update path length this->_path_length += h; - this->_s += h; + this->_path_length_per_surface += h; } template ::state:: advance_jacobian(const detray::stepping::config& cfg) { + + // Immediately exit if it is for simulation + if (cfg.is_simulation) { + return; + } + /// The calculations are based on ATL-SOFT-PUB-2009-002. The update of the /// Jacobian matrix is requires only the calculation of eq. 17 and 18. /// Since the terms of eq. 18 are currently 0, this matrix is not needed @@ -414,22 +444,29 @@ detray::rk_stepper class array_t> -DETRAY_HOST_DEVICE void -detray::rk_stepper::state:: - add_multiple_scattering_covariance( - const detray::stepping::config& cfg) {} +DETRAY_HOST_DEVICE void detray::rk_stepper< + magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t, + inspector_t, array_t>::state::add_multiple_scattering_covariance() { + + // Implement the thick scatterer method of Eq (4.103 - 111) of "Pattern + // Recognition, Tracking and Vertex Reconstruction in Particle Detectors" + + // Variance per unit length of the projected scattering angle + // interaction_type(). + + // 4X4 Identity matrix + const matrix_type<4, 4> E = matrix_operator().template zero<4, 4>(); +} template class array_t> -DETRAY_HOST_DEVICE auto -detray::rk_stepper::state:: - evaluate_dqopds(const std::size_t i, - const typename transform3_t::scalar_type h, - const scalar dqopds_prev, - const detray::stepping::config& cfg) -> +DETRAY_HOST_DEVICE auto detray::rk_stepper< + magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t, + inspector_t, + array_t>::state::evaluate_dqopds(const std::size_t i, + const typename transform3_t::scalar_type h, + const scalar dqopds_prev) -> typename transform3_t::scalar_type { const auto& track = this->_track; @@ -441,20 +478,21 @@ detray::rk_stepperdqopds(sd.qop[i]); } } @@ -645,7 +683,7 @@ DETRAY_HOST_DEVICE bool detray::rk_stepper< // qop should be recalcuated at every point // Reference: Eq (84) of https://doi.org/10.1016/0029-554X(81)90063-X - sd.dqopds[0u] = stepping.evaluate_dqopds(0u, 0.f, 0.f, cfg); + sd.dqopds[0u] = stepping.evaluate_dqopds(0u, 0.f, 0.f); sd.dtds[0u] = stepping.evaluate_dtds(sd.b_first, 0u, 0.f, vector3{0.f, 0.f, 0.f}, sd.qop[0u]); @@ -665,7 +703,7 @@ DETRAY_HOST_DEVICE bool detray::rk_stepper< sd.b_middle[2] = bvec1[2]; sd.dqopds[1u] = - stepping.evaluate_dqopds(1u, half_h, sd.dqopds[0u], cfg); + stepping.evaluate_dqopds(1u, half_h, sd.dqopds[0u]); sd.dtds[1u] = stepping.evaluate_dtds(sd.b_middle, 1u, half_h, sd.dtds[0u], sd.qop[1u]); @@ -673,7 +711,7 @@ DETRAY_HOST_DEVICE bool detray::rk_stepper< // qop should be recalcuated at every point // Reference: Eq (84) of https://doi.org/10.1016/0029-554X(81)90063-X sd.dqopds[2u] = - stepping.evaluate_dqopds(2u, half_h, sd.dqopds[1u], cfg); + stepping.evaluate_dqopds(2u, half_h, sd.dqopds[1u]); sd.dtds[2u] = stepping.evaluate_dtds(sd.b_middle, 2u, half_h, sd.dtds[1u], sd.qop[2u]); @@ -686,7 +724,7 @@ DETRAY_HOST_DEVICE bool detray::rk_stepper< sd.b_last[1] = bvec2[1]; sd.b_last[2] = bvec2[2]; - sd.dqopds[3u] = stepping.evaluate_dqopds(3u, h, sd.dqopds[2u], cfg); + sd.dqopds[3u] = stepping.evaluate_dqopds(3u, h, sd.dqopds[2u]); sd.dtds[3u] = stepping.evaluate_dtds(sd.b_last, 3u, h, sd.dtds[2u], sd.qop[3u]); diff --git a/core/include/detray/propagator/stepping_config.hpp b/core/include/detray/propagator/stepping_config.hpp index 48be594b1..4a1d03997 100644 --- a/core/include/detray/propagator/stepping_config.hpp +++ b/core/include/detray/propagator/stepping_config.hpp @@ -35,15 +35,19 @@ struct config { scalar_t path_limit{5.f * unit::m}; /// Maximum number of Runge-Kutta step trials std::size_t max_rk_updates{10000u}; - /// Include multiple scattering -> only for simulation + /* + /// Scatter trajectory for simulation bool do_scattering{false}; /// Use mean energy loss (Bethe) /// if false, most probable energy loss (Landau) will be used bool use_mean_loss{true}; + */ /// Use eloss gradient in error propagation bool use_eloss_gradient{false}; /// Use b field gradient in error propagation bool use_field_gradient{false}; + /// Is this for reconstruction or simulatoin + bool is_simulation{false}; }; } // namespace detray::stepping diff --git a/tests/integration_tests/cpu/propagator/propagator.cpp b/tests/integration_tests/cpu/propagator/propagator.cpp index 43bf28e15..0e0c1f871 100644 --- a/tests/integration_tests/cpu/propagator/propagator.cpp +++ b/tests/integration_tests/cpu/propagator/propagator.cpp @@ -71,7 +71,8 @@ struct helix_inspector : actor { } // Nothing has happened yet (first call of actor chain) - if (stepping.path_length() < tol || stepping._s < tol) { + if (stepping.path_length() < tol || + stepping._path_length_per_surface < tol) { return; } @@ -94,21 +95,21 @@ struct helix_inspector : actor { detail::helix hlx(free_params, &b); - const auto true_pos = hlx(stepping._s); + const auto true_pos = hlx(stepping._path_length_per_surface); - const point3 relative_error{1.f / stepping._s * + const point3 relative_error{1.f / stepping._path_length_per_surface * (stepping().pos() - true_pos)}; ASSERT_NEAR(getter::norm(relative_error), 0.f, tol); - auto true_J = hlx.jacobian(stepping._s); + auto true_J = hlx.jacobian(stepping._path_length_per_surface); for (unsigned int i = 0u; i < e_free_size; i++) { for (unsigned int j = 0u; j < e_free_size; j++) { ASSERT_NEAR( matrix_operator().element(stepping._jac_transport, i, j), matrix_operator().element(true_J, i, j), - stepping._s * tol * 10.f); + stepping._path_length_per_surface * tol * 10.f); } } } From d34fc4b99d7ec3c90286eee9aaa8cbd170d37e47 Mon Sep 17 00:00:00 2001 From: beomki-yeo Date: Sat, 17 Feb 2024 22:16:57 +0100 Subject: [PATCH 08/13] Backup --- .../propagator/actors/parameter_resetter.hpp | 6 ++ .../actors/parameter_transporter.hpp | 9 +- .../detray/propagator/base_stepper.hpp | 11 +++ core/include/detray/propagator/rk_stepper.hpp | 3 +- core/include/detray/propagator/rk_stepper.ipp | 95 +++++++++++++++---- .../cpu/propagator/jacobian_validation.cpp | 20 ++-- 6 files changed, 116 insertions(+), 28 deletions(-) diff --git a/core/include/detray/propagator/actors/parameter_resetter.hpp b/core/include/detray/propagator/actors/parameter_resetter.hpp index b9214ffc9..f33a3e837 100644 --- a/core/include/detray/propagator/actors/parameter_resetter.hpp +++ b/core/include/detray/propagator/actors/parameter_resetter.hpp @@ -51,12 +51,18 @@ struct parameter_resetter : actor { // Reset the path length stepping._path_length_per_surface = 0; + // Reset the qop at the surface + stepping._qop_i = stepping().qop(); + // Reset jacobian coordinate transformation at the current surface stepping._jac_to_global = local_coordinate.bound_to_free_jacobian( trf3, mask, stepping._bound_params.vector()); // Reset jacobian transport to identity matrix matrix_operator().set_identity(stepping._jac_transport); + + // Reset the joint covariance + matrix_operator().set_zero(stepping._joint_cov); } }; diff --git a/core/include/detray/propagator/actors/parameter_transporter.hpp b/core/include/detray/propagator/actors/parameter_transporter.hpp index 91921b946..e177f6f5d 100644 --- a/core/include/detray/propagator/actors/parameter_transporter.hpp +++ b/core/include/detray/propagator/actors/parameter_transporter.hpp @@ -12,6 +12,7 @@ #include "detray/definitions/track_parametrization.hpp" #include "detray/geometry/surface.hpp" #include "detray/propagator/base_actor.hpp" +#include "detray/propagator/propagation_config.hpp" namespace detray { @@ -119,11 +120,9 @@ struct parameter_transporter : actor { matrix_operator().transpose(stepping._full_jacobian); } - /* - // Update covariance with the multiple scattering - stepping.add_multiple_scattering_covariance(); - */ - + // Add the multiple scattering term + new_cov = new_cov + stepping._joint_cov; + // Calculate surface-to-surface covariance transport stepping._bound_params.set_covariance(new_cov); } diff --git a/core/include/detray/propagator/base_stepper.hpp b/core/include/detray/propagator/base_stepper.hpp index e9f9a364a..3672344f3 100644 --- a/core/include/detray/propagator/base_stepper.hpp +++ b/core/include/detray/propagator/base_stepper.hpp @@ -76,6 +76,7 @@ class base_stepper { struct state { /// Sets track parameters. + /// @TODO: Remove the free track parameter input if possible DETRAY_HOST_DEVICE state(const free_track_parameters_type &t) : _track(t) {} @@ -93,9 +94,12 @@ class base_stepper { sf.template visit_mask< typename parameter_resetter::kernel>( sf.transform(ctx), *this); + + _qop_i = bound_params.qop(); } /// free track parameter + /// @TODO: Remove the covariance if possible free_track_parameters_type _track; /// Full jacobian @@ -113,6 +117,10 @@ class base_stepper { /// bound covariance bound_track_parameters_type _bound_params; + /// joint covariance (for multiple scattering) + bound_matrix _joint_cov = + matrix_operator().template zero(); + /// @returns track parameters - const access DETRAY_HOST_DEVICE free_track_parameters_type &operator()() { return _track; } @@ -148,6 +156,9 @@ class base_stepper { /// The particle pdg int _pdg = 13; // default muon + /// qop at the initial surface + scalar_type _qop_i{0.f}; + /// Random device random_device_t rand_device; diff --git a/core/include/detray/propagator/rk_stepper.hpp b/core/include/detray/propagator/rk_stepper.hpp index 890cf20a1..a68e38090 100644 --- a/core/include/detray/propagator/rk_stepper.hpp +++ b/core/include/detray/propagator/rk_stepper.hpp @@ -107,7 +107,8 @@ class rk_stepper final /// Update the covariance from multiple scattering DETRAY_HOST_DEVICE - inline void add_multiple_scattering_covariance(); + inline void add_multiple_scattering_covariance( + const detray::stepping::config& cfg); /// evaulate dqopds for a given step size and material DETRAY_HOST_DEVICE diff --git a/core/include/detray/propagator/rk_stepper.ipp b/core/include/detray/propagator/rk_stepper.ipp index 774c20e3c..f4d604d26 100644 --- a/core/include/detray/propagator/rk_stepper.ipp +++ b/core/include/detray/propagator/rk_stepper.ipp @@ -70,18 +70,16 @@ detray::rk_stepper_path_length_per_surface, mat, this->_pdg, - this->_mass, qop, track.charge()); + h, mat, this->_pdg, this->_mass, qop, track.charge()); const scalar_type e_loss_sigma = interaction_type().compute_energy_loss_landau_sigma( - this->_path_length_per_surface, mat, this->_pdg, - this->_mass, qop, track.charge()); + h, mat, this->_pdg, this->_mass, qop, track.charge()); // Get the new momentum - const auto new_mom = attenuate(e_loss_mpv, e_loss_sigma, - this->mass, track.charge() / qop, - this->rand_device.generator); + const auto new_mom = random_scatterer().attenuate( + e_loss_mpv, e_loss_sigma, this->mass, track.charge() / qop, + this->rand_device.generator); qop = track.charge() / new_mom; } @@ -444,18 +442,85 @@ detray::rk_stepper class array_t> -DETRAY_HOST_DEVICE void detray::rk_stepper< - magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t, - inspector_t, array_t>::state::add_multiple_scattering_covariance() { +DETRAY_HOST_DEVICE void +detray::rk_stepper::state:: + add_multiple_scattering_covariance( + const detray::stepping::config& cfg) { + + // Return if there is no material + if (this->mat == vacuum()) { + return; + } + + // Return if it is for simulation + if (cfg.is_simulation) { + return; + } // Implement the thick scatterer method of Eq (4.103 - 111) of "Pattern // Recognition, Tracking and Vertex Reconstruction in Particle Detectors" // Variance per unit length of the projected scattering angle - // interaction_type(). + const scalar s1 = interaction_type().compute_multiple_scattering_theta0( + 1.f / this->mat.X0(), this->_pdg, this->_mass, this->_qop_i, + this->_track.charge()); + const scalar s2 = interaction_type().compute_multiple_scattering_theta0( + 1.f / this->mat.X0(), this->_pdg, this->_mass, this->_track.qop(), + this->_track.charge()); + const scalar variance = s1 * s2; // 4X4 Identity matrix - const matrix_type<4, 4> E = matrix_operator().template zero<4, 4>(); + matrix_type<4, 4> E = matrix_operator().template zero<4, 4>(); + + const scalar L = this->_path_length_per_surface; + const scalar half_L2 = L * L / 2.f; + const scalar third_L3 = L * L * L / 3.f; + + // Set E (Eq. 4.107). Note that we flip the index of (0,1) <-> (2,3) + getter::element(E, 0, 0) = third_L3; + getter::element(E, 1, 1) = third_L3; + getter::element(E, 2, 2) = L; + getter::element(E, 3, 3) = L; + getter::element(E, 0, 2) = half_L2; + getter::element(E, 1, 3) = half_L2; + getter::element(E, 2, 0) = half_L2; + getter::element(E, 3, 1) = half_L2; + + const matrix_type<3, 3> C2G = + mat_helper().curvilinear_to_global(this->_track.dir()); + // Take the bottom-right 2x2 block + const matrix_type<2, 2> block = + matrix_operator().template block<2, 2>(C2G, 0, 0); + + // Jacobian for (Eq. 4.108). Note that we flip the index of (0,1) <-> (2,3) + matrix_type<4, 4> T = matrix_operator().template identity<4, 4>(); + matrix_operator().template set_block<2, 2>(T, block, 2, 2); + + E = T * E * matrix_operator().transpose(T); + + // Jacobian d(phi,theta)/d(tx,ty). Note that it is different from + // (Eq. 4.109) due to the different coordinate system + matrix_type<4, 4> J = matrix_operator().template identity<4, 4>(); + + const auto t = this->_track.dir(); + const scalar_type phi = getter::phi(t); + const scalar_type theta = getter::theta(t); + const scalar_type sin_phi = math::sin(phi); + const scalar_type cos_phi = math::cos(phi); + const scalar_type sin_theta = math::sin(theta); + const scalar_type cos_theta = math::cos(theta); + + // -sin(phi)/sin(theta) + getter::element(J, 2, 2) = -sin_phi / sin_theta; + getter::element(J, 2, 3) = -cos_phi / sin_theta; + getter::element(J, 3, 2) = -cos_phi * cos_theta; + getter::element(J, 3, 3) = -sin_phi * cos_theta; + + E = J * E * matrix_operator().transpose(J); + + // Set the joint covariance + matrix_operator().set_block<4, 4>(this->_joint_cov, E, 0, 0); } template ::state evaluate_bound_param( const typename propagator_t::detector_type& det, const field_t& field, const scalar overstep_tolerance, const scalar on_surface_tolerance, const scalar rk_tolerance, const scalar constraint_step, - bool use_field_gradient, bool do_inspect = false) { + bool use_field_gradient, bool do_scatter, bool do_inspect = false) { // Propagator is built from the stepper and navigator propagation::config cfg{}; @@ -354,6 +354,7 @@ bound_getter::state evaluate_bound_param( cfg.stepping.rk_error_tol = rk_tolerance; cfg.stepping.use_eloss_gradient = true; cfg.stepping.use_field_gradient = use_field_gradient; + cfg.stepping.is_simulation = do_scatter; propagator_t p(cfg); // Actor states @@ -625,10 +626,13 @@ void evaluate_jacobian_difference( det.volumes()[0u].set_material(volume_mat); + const bool do_scatter = true; + + // No scattering for jacobian comparsion study auto bound_getter = evaluate_bound_param( detector_length, track, det, field, overstep_tolerance, on_surface_tolerance, rk_tolerance, constraint_step, use_field_gradient, - do_inspect); + !do_scatter, do_inspect); const auto reference_param = bound_getter.m_param_departure; const auto final_param = bound_getter.m_param_destination; @@ -754,10 +758,13 @@ void evaluate_covariance_transport( track_copy.set_covariance(ini_cov); + const bool do_scatter = true; + + // Reference track (Without Scattering!) auto bound_getter = evaluate_bound_param( detector_length, track_copy, det, field, overstep_tolerance, - on_surface_tolerance, rk_tolerance, constraint_step, - use_field_gradient); + on_surface_tolerance, rk_tolerance, constraint_step, use_field_gradient, + !do_scatter); const auto reference_param = bound_getter.m_param_departure; const auto ini_vec = reference_param.vector(); @@ -793,10 +800,11 @@ void evaluate_covariance_transport( auto smeared_track = track_copy; smeared_track.set_vector(smeared_ini_vec); + // Smeared track (With scattering!) auto smeared_bound_getter = evaluate_bound_param( detector_length, smeared_track, det, field, overstep_tolerance, - on_surface_tolerance, rk_tolerance, constraint_step, - use_field_gradient); + on_surface_tolerance, rk_tolerance, constraint_step, use_field_gradient, + do_scatter); // Get smeared final bound vector bound_vector_type smeared_fin_vec = From 95fc25cfcfb412e4127fe70b602591eac2b59794 Mon Sep 17 00:00:00 2001 From: beomki-yeo Date: Sat, 17 Feb 2024 22:48:51 +0100 Subject: [PATCH 09/13] Backup --- .../detray/propagator/line_stepper.hpp | 2 +- core/include/detray/propagator/rk_stepper.hpp | 4 ++-- core/include/detray/propagator/rk_stepper.ipp | 23 +++++++++++-------- .../detray/propagator/stepping_config.hpp | 7 ++++-- .../cpu/propagator/jacobian_validation.cpp | 22 ++++++++++++------ 5 files changed, 36 insertions(+), 22 deletions(-) diff --git a/core/include/detray/propagator/line_stepper.hpp b/core/include/detray/propagator/line_stepper.hpp index 1b0581cb2..0d6e530a0 100644 --- a/core/include/detray/propagator/line_stepper.hpp +++ b/core/include/detray/propagator/line_stepper.hpp @@ -81,7 +81,7 @@ class line_stepper final } DETRAY_HOST_DEVICE - inline void add_multiple_scattering_covariance() { return; } + inline void calculate_ms_covariance() { return; } DETRAY_HOST_DEVICE inline vector3 dtds() const { return vector3{0.f, 0.f, 0.f}; } diff --git a/core/include/detray/propagator/rk_stepper.hpp b/core/include/detray/propagator/rk_stepper.hpp index a68e38090..9601061fa 100644 --- a/core/include/detray/propagator/rk_stepper.hpp +++ b/core/include/detray/propagator/rk_stepper.hpp @@ -107,8 +107,8 @@ class rk_stepper final /// Update the covariance from multiple scattering DETRAY_HOST_DEVICE - inline void add_multiple_scattering_covariance( - const detray::stepping::config& cfg); + inline void calculate_ms_covariance( + const detray::stepping::config& cfg = {}); /// evaulate dqopds for a given step size and material DETRAY_HOST_DEVICE diff --git a/core/include/detray/propagator/rk_stepper.ipp b/core/include/detray/propagator/rk_stepper.ipp index f4d604d26..839940575 100644 --- a/core/include/detray/propagator/rk_stepper.ipp +++ b/core/include/detray/propagator/rk_stepper.ipp @@ -46,7 +46,7 @@ detray::rk_stepper::value) { // Scatter if it is for the simulation - if (cfg.is_simulation) { + if (cfg.do_scatter) { const scalar_type projected_scattering_angle = interaction_type().compute_multiple_scattering_theta0( h / mat.X0(), this->_pdg, this->_mass, qop, @@ -59,7 +59,7 @@ detray::rk_stepper& cfg) { // Immediately exit if it is for simulation - if (cfg.is_simulation) { + if (cfg.do_scatter || !cfg.do_covariance_transport) { return; } @@ -445,16 +445,15 @@ template ::state:: - add_multiple_scattering_covariance( - const detray::stepping::config& cfg) { + calculate_ms_covariance(const detray::stepping::config& cfg) { // Return if there is no material - if (this->mat == vacuum()) { + if (this->_mat == vacuum()) { return; } // Return if it is for simulation - if (cfg.is_simulation) { + if (cfg.do_scatter || !cfg.do_covariance_transport) { return; } @@ -463,10 +462,10 @@ detray::rk_steppermat.X0(), this->_pdg, this->_mass, this->_qop_i, + 1.f / this->_mat.X0(), this->_pdg, this->_mass, this->_qop_i, this->_track.charge()); const scalar s2 = interaction_type().compute_multiple_scattering_theta0( - 1.f / this->mat.X0(), this->_pdg, this->_mass, this->_track.qop(), + 1.f / this->_mat.X0(), this->_pdg, this->_mass, this->_track.qop(), this->_track.charge()); const scalar variance = s1 * s2; @@ -486,6 +485,7 @@ detray::rk_stepper C2G = mat_helper().curvilinear_to_global(this->_track.dir()); @@ -520,7 +520,7 @@ detray::rk_stepper(this->_joint_cov, E, 0, 0); + matrix_operator().template set_block<4, 4>(this->_joint_cov, E, 0, 0); } template ::state evaluate_bound_param( const typename propagator_t::detector_type& det, const field_t& field, const scalar overstep_tolerance, const scalar on_surface_tolerance, const scalar rk_tolerance, const scalar constraint_step, - bool use_field_gradient, bool do_scatter, bool do_inspect = false) { + bool use_field_gradient, bool do_scatter, bool do_covariance_transport, + bool do_inspect = false) { // Propagator is built from the stepper and navigator propagation::config cfg{}; @@ -354,7 +355,8 @@ bound_getter::state evaluate_bound_param( cfg.stepping.rk_error_tol = rk_tolerance; cfg.stepping.use_eloss_gradient = true; cfg.stepping.use_field_gradient = use_field_gradient; - cfg.stepping.is_simulation = do_scatter; + cfg.stepping.do_scatter = do_scatter; + cfg.stepping.do_covariance_transport = do_covariance_transport; propagator_t p(cfg); // Actor states @@ -396,6 +398,9 @@ get_displaced_bound_vector( cfg.navigation.overstep_tolerance = overstep_tolerance; cfg.navigation.on_surface_tolerance = on_surface_tolerance; cfg.stepping.rk_error_tol = rk_tolerance; + // we don't have to scatering or covariance transport + cfg.stepping.do_scatter = false; + cfg.stepping.do_covariance_transport = false; // Propagator is built from the stepper and navigator propagator_t p(cfg); @@ -627,12 +632,14 @@ void evaluate_jacobian_difference( det.volumes()[0u].set_material(volume_mat); const bool do_scatter = true; + const bool do_covariance_transport = true; // No scattering for jacobian comparsion study + // But perform covariance transport auto bound_getter = evaluate_bound_param( detector_length, track, det, field, overstep_tolerance, on_surface_tolerance, rk_tolerance, constraint_step, use_field_gradient, - !do_scatter, do_inspect); + !do_scatter, do_covariance_transport, do_inspect); const auto reference_param = bound_getter.m_param_departure; const auto final_param = bound_getter.m_param_destination; @@ -759,12 +766,13 @@ void evaluate_covariance_transport( track_copy.set_covariance(ini_cov); const bool do_scatter = true; + const bool do_covariance_transport = true; - // Reference track (Without Scattering!) + // Reference track (Without Scattering but with covariance transport) auto bound_getter = evaluate_bound_param( detector_length, track_copy, det, field, overstep_tolerance, on_surface_tolerance, rk_tolerance, constraint_step, use_field_gradient, - !do_scatter); + !do_scatter, do_covariance_transport); const auto reference_param = bound_getter.m_param_departure; const auto ini_vec = reference_param.vector(); @@ -800,11 +808,11 @@ void evaluate_covariance_transport( auto smeared_track = track_copy; smeared_track.set_vector(smeared_ini_vec); - // Smeared track (With scattering!) + // Smeared track (With scattering and no covariance transport) auto smeared_bound_getter = evaluate_bound_param( detector_length, smeared_track, det, field, overstep_tolerance, on_surface_tolerance, rk_tolerance, constraint_step, use_field_gradient, - do_scatter); + do_scatter, !do_covariance_transport); // Get smeared final bound vector bound_vector_type smeared_fin_vec = From 5f0425bb5c02ef874684c34e5b818a0ccb7b9548 Mon Sep 17 00:00:00 2001 From: beomki-yeo Date: Sat, 17 Feb 2024 23:08:53 +0100 Subject: [PATCH 10/13] Backup --- tests/scripts/run_jacobian_validation.sh | 80 +++++++++++++++--------- 1 file changed, 49 insertions(+), 31 deletions(-) diff --git a/tests/scripts/run_jacobian_validation.sh b/tests/scripts/run_jacobian_validation.sh index ed362cef8..6eb66d7b6 100644 --- a/tests/scripts/run_jacobian_validation.sh +++ b/tests/scripts/run_jacobian_validation.sh @@ -27,13 +27,16 @@ log10_on_surface_tol=-3 # Monte-Carlo seed mc_seed=0 +# Skip the first phase +skip_first_phase=false + # Skip the second phase skip_second_phase=false # Verbose level verbose_level=1 -while getopts "hd:n:t:p:q:c:i:s:r:v:" arg; do +while getopts "hd:n:t:p:q:c:i:s:f:r:v:" arg; do case $arg in h) echo "" @@ -48,6 +51,7 @@ while getopts "hd:n:t:p:q:c:i:s:r:v:" arg; do echo "-c " echo "-i " echo "-s " + echo "-f " echo "-r " echo "-v " echo "" @@ -86,6 +90,10 @@ while getopts "hd:n:t:p:q:c:i:s:r:v:" arg; do mc_seed=$OPTARG echo "Monte-Carlo seed: ${mc_seed}" ;; + f) + skip_first_phase=$OPTARG + echo "Skip the first phase: ${skip_first_phase}" + ;; r) skip_second_phase=$OPTARG echo "Skip the second phase: ${skip_second_phase}" @@ -108,35 +116,39 @@ fi # RK tolerance iteration # ########################## -echo "Starting rk toleracne iteration..." - -# Remove the old directories -for (( i=0; i < ${n_threads}; ++i )) -do - rm -rf ${PWD}/thread_${i} -done -rm -rf ${PWD}/merged - -for (( i=0; i < ${n_threads}; ++i )) -do - n_skips=`expr ${i} \* ${n_tracks_per_thread}` +if [ "$skip_first_phase" = false ] ; then - command_rk_tolerance="${dir}/detray_integration_test_jacobian_validation \ - --output-directory=thread_${i} \ - --rk-tolerance-iterate-mode=true \ - --n-tracks=${n_tracks_per_thread} \ - --n-skips=${n_skips} \ - --log10-min-rk-tolerance=${log10_min_rk_tol} \ - --log10-max-rk-tolerance=${log10_max_rk_tol} \ - --log10-helix-tolerance=${log10_helix_tol} \ - --log10-on-surface-tolerance=${log10_on_surface_tol} \ - --mc-seed=${mc_seed} \ - --verbose-level=${verbose_level}" - ${command_rk_tolerance} & -done -wait - -echo "Finished rk toleracne iteration" + echo "Starting rk toleracne iteration..." + + # Remove the old directories + for (( i=0; i < ${n_threads}; ++i )) + do + rm -rf ${PWD}/thread_${i} + done + rm -rf ${PWD}/merged + + for (( i=0; i < ${n_threads}; ++i )) + do + n_skips=`expr ${i} \* ${n_tracks_per_thread}` + + command_rk_tolerance="${dir}/detray_integration_test_jacobian_validation \ + --output-directory=thread_${i} \ + --rk-tolerance-iterate-mode=true \ + --n-tracks=${n_tracks_per_thread} \ + --n-skips=${n_skips} \ + --log10-min-rk-tolerance=${log10_min_rk_tol} \ + --log10-max-rk-tolerance=${log10_max_rk_tol} \ + --log10-helix-tolerance=${log10_helix_tol} \ + --log10-on-surface-tolerance=${log10_on_surface_tol} \ + --mc-seed=${mc_seed} \ + --verbose-level=${verbose_level}" + ${command_rk_tolerance} & + done + wait + + echo "Finished rk toleracne iteration" + +fi ##################################### # Jacobi validation & Cov transport # @@ -209,7 +221,7 @@ echo "Finished merging Csv files" cd ${output_dir} -if [ "$skip_second_phase" = false ] ; then +if [ "$skip_first_phase" = false ] && [ "$skip_second_phase" = false ]; then # Run rk_tolerance_comparision.C root -q '../../../tests/validation/root/rk_tolerance_comparison.C+O('${log10_min_rk_tol}','${log10_max_rk_tol}')' @@ -219,7 +231,13 @@ if [ "$skip_second_phase" = false ] ; then # Run covariance_validation.C root -q -l ../../../tests/validation/root/covariance_validation.C+O -else + + elif [ "$skip_first_phase" = true ] && [ "$skip_second_phase" = false ]; then + + # Run covariance_validation.C + root -q -l ../../../tests/validation/root/covariance_validation.C+O + + elif [ "$skip_first_phase" = false ] && [ "$skip_second_phase" = true ]; then # Run rk_tolerance_comparision.C root '../../../tests/validation/root/rk_tolerance_comparison.C+O('${log10_min_rk_tol}','${log10_max_rk_tol}')' From ff42a9ad6ce4a020178768ab6e01bc3dc3f84c0d Mon Sep 17 00:00:00 2001 From: beomki-yeo Date: Sun, 18 Feb 2024 00:33:06 +0100 Subject: [PATCH 11/13] Backup --- .../actors/parameter_transporter.hpp | 37 ++++++++++++++++--- .../detray/propagator/line_stepper.hpp | 6 ++- core/include/detray/propagator/rk_stepper.hpp | 6 ++- core/include/detray/propagator/rk_stepper.ipp | 28 ++++++++------ .../cpu/propagator/jacobian_validation.cpp | 11 ++++-- 5 files changed, 63 insertions(+), 25 deletions(-) diff --git a/core/include/detray/propagator/actors/parameter_transporter.hpp b/core/include/detray/propagator/actors/parameter_transporter.hpp index e177f6f5d..bc703c63a 100644 --- a/core/include/detray/propagator/actors/parameter_transporter.hpp +++ b/core/include/detray/propagator/actors/parameter_transporter.hpp @@ -14,12 +14,18 @@ #include "detray/propagator/base_actor.hpp" #include "detray/propagator/propagation_config.hpp" +#include + namespace detray { template struct parameter_transporter : actor { - struct state {}; + struct state { + // @TODO remove once we can call stepping config... + bool do_scatter = false; + bool do_covariance_transport = true; + }; /// Mask store visitor struct kernel { @@ -59,7 +65,8 @@ struct parameter_transporter : actor { typename propagator_state_t> DETRAY_HOST_DEVICE inline void operator()( const mask_group_t& mask_group, const index_t& index, - const transform3_type& trf3, propagator_state_t& propagation) { + const transform3_type& trf3, propagator_state_t& propagation, + state& actor_state) { // Stepper and Navigator states auto& stepping = propagation._stepping; @@ -120,16 +127,33 @@ struct parameter_transporter : actor { matrix_operator().transpose(stepping._full_jacobian); } - // Add the multiple scattering term - new_cov = new_cov + stepping._joint_cov; + // Calculate multiple scattering term + // @TODO: Take stepping::config instead + const auto joint_cov = stepping.calculate_ms_covariance( + actor_state.do_scatter, actor_state.do_covariance_transport); + // Add the multiple scattering term + // new_cov = new_cov + stepping._joint_cov; + new_cov = new_cov + joint_cov; + + /* + std::cout << "Joint Cov" << std::endl; + for (int i = 0; i < 6; i++) { + for (int j = 0; j < 6; j++) { + std::cout << getter::element(new_cov, i, j) << " "; + } + std::cout << std::endl; + } + std::cout << std::endl; + */ + // Calculate surface-to-surface covariance transport stepping._bound_params.set_covariance(new_cov); } }; template - DETRAY_HOST_DEVICE void operator()(state& /*actor_state*/, + DETRAY_HOST_DEVICE void operator()(state& actor_state, propagator_state_t& propagation) const { const auto& navigation = propagation._navigation; @@ -145,7 +169,8 @@ struct parameter_transporter : actor { // Surface const auto sf = navigation.get_surface(); - sf.template visit_mask(sf.transform(ctx), propagation); + sf.template visit_mask(sf.transform(ctx), propagation, + actor_state); // Set surface link propagation._stepping._bound_params.set_surface_link(sf.barcode()); diff --git a/core/include/detray/propagator/line_stepper.hpp b/core/include/detray/propagator/line_stepper.hpp index 0d6e530a0..497eed101 100644 --- a/core/include/detray/propagator/line_stepper.hpp +++ b/core/include/detray/propagator/line_stepper.hpp @@ -32,6 +32,7 @@ class line_stepper final typename base_type::free_track_parameters_type; using bound_track_parameters_type = typename base_type::bound_track_parameters_type; + using bound_matrix = typename bound_track_parameters_type::covariance_type; using matrix_operator = typename base_type::matrix_operator; using size_type = typename matrix_operator::size_ty; using vector3 = typename line_stepper::transform3_type::vector3; @@ -81,7 +82,10 @@ class line_stepper final } DETRAY_HOST_DEVICE - inline void calculate_ms_covariance() { return; } + inline bound_matrix calculate_ms_covariance( + bool /*do_scatter*/, bool /*do_covariance_transport*/) const { + return matrix_operator().template zero<6, 6>(); + } DETRAY_HOST_DEVICE inline vector3 dtds() const { return vector3{0.f, 0.f, 0.f}; } diff --git a/core/include/detray/propagator/rk_stepper.hpp b/core/include/detray/propagator/rk_stepper.hpp index 9601061fa..dbaa2c8cc 100644 --- a/core/include/detray/propagator/rk_stepper.hpp +++ b/core/include/detray/propagator/rk_stepper.hpp @@ -49,6 +49,7 @@ class rk_stepper final typename base_type::free_track_parameters_type; using bound_track_parameters_type = typename base_type::bound_track_parameters_type; + using bound_matrix = typename bound_track_parameters_type::covariance_type; using magnetic_field_type = magnetic_field_t; using size_type = typename matrix_operator::size_ty; template @@ -106,9 +107,10 @@ class rk_stepper final const stepping::config& cfg = {}); /// Update the covariance from multiple scattering + /// @TODO Take stepping::config DETRAY_HOST_DEVICE - inline void calculate_ms_covariance( - const detray::stepping::config& cfg = {}); + inline bound_matrix calculate_ms_covariance( + bool do_scatter, bool do_covariance_transport) const; /// evaulate dqopds for a given step size and material DETRAY_HOST_DEVICE diff --git a/core/include/detray/propagator/rk_stepper.ipp b/core/include/detray/propagator/rk_stepper.ipp index 839940575..110526fa0 100644 --- a/core/include/detray/propagator/rk_stepper.ipp +++ b/core/include/detray/propagator/rk_stepper.ipp @@ -78,7 +78,7 @@ detray::rk_stepper().attenuate( - e_loss_mpv, e_loss_sigma, this->mass, track.charge() / qop, + e_loss_mpv, e_loss_sigma, this->_mass, track.charge() / qop, this->rand_device.generator); qop = track.charge() / new_mom; @@ -442,19 +442,24 @@ detray::rk_stepper class array_t> -DETRAY_HOST_DEVICE void -detray::rk_stepper::state:: - calculate_ms_covariance(const detray::stepping::config& cfg) { +DETRAY_HOST_DEVICE auto detray::rk_stepper< + magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t, + inspector_t, + array_t>::state::calculate_ms_covariance(bool do_scatter, + bool do_covariance_transport) + const -> bound_matrix { + + // Joint covariance + bound_matrix joint_cov = matrix_operator().template zero<6,6>(); // Return if there is no material if (this->_mat == vacuum()) { - return; + return joint_cov; } // Return if it is for simulation - if (cfg.do_scatter || !cfg.do_covariance_transport) { - return; + if (do_scatter || do_covariance_transport) { + return joint_cov; } // Implement the thick scatterer method of Eq (4.103 - 111) of "Pattern @@ -520,7 +525,9 @@ detray::rk_stepper(this->_joint_cov, E, 0, 0); + matrix_operator().template set_block<4, 4>(joint_cov, E, 0, 0); + + return joint_cov; } template ::state evaluate_bound_param( propagator_t p(cfg); // Actor states - parameter_transporter::state transporter_state{}; + parameter_transporter::state transporter_state{ + cfg.stepping.do_scatter, cfg.stepping.do_covariance_transport}; bound_getter::state bound_getter_state{}; bound_getter_state.m_min_path_length = detector_length * 0.75f; parameter_resetter::state resetter_state{}; @@ -414,7 +416,8 @@ get_displaced_bound_vector( typename propagator_t::state dstate(dparam, field, det); // Actor states - parameter_transporter::state transporter_state{}; + parameter_transporter::state transporter_state{ + cfg.stepping.do_scatter, cfg.stepping.do_covariance_transport}; parameter_resetter::state resetter_state{}; bound_getter::state bound_getter_state{}; bound_getter_state.m_min_path_length = detector_length * 0.75f; @@ -1522,10 +1525,10 @@ int main(int argc, char** argv) { // Stepper types using const_field_stepper_t = rk_stepper, - stepper_default_policy>; + stepper_default_policy, stepping::default_random_device>; using inhom_field_stepper_t = rk_stepper, - stepper_default_policy>; + stepper_default_policy, stepping::default_random_device>; // Make four propagators for each case using const_field_rect_propagator_t = From d3181a443088ca3850c4d996d74b4bc5766b8396 Mon Sep 17 00:00:00 2001 From: beomki-yeo Date: Sun, 18 Feb 2024 01:11:45 +0100 Subject: [PATCH 12/13] Backup --- .gitignore | 2 +- .../actors/parameter_transporter.hpp | 7 +- .../cpu/propagator/jacobian_validation.cpp | 100 ++++++++++-------- tests/scripts/run_jacobian_validation.sh | 16 ++- 4 files changed, 76 insertions(+), 49 deletions(-) diff --git a/.gitignore b/.gitignore index 4e3d751dc..0d3e2fbc3 100644 --- a/.gitignore +++ b/.gitignore @@ -15,4 +15,4 @@ **/*.d **/*.so **/*C_ACLiC_dict* -**/therad_* \ No newline at end of file +**/thread_* \ No newline at end of file diff --git a/core/include/detray/propagator/actors/parameter_transporter.hpp b/core/include/detray/propagator/actors/parameter_transporter.hpp index bc703c63a..a2340eee2 100644 --- a/core/include/detray/propagator/actors/parameter_transporter.hpp +++ b/core/include/detray/propagator/actors/parameter_transporter.hpp @@ -127,6 +127,7 @@ struct parameter_transporter : actor { matrix_operator().transpose(stepping._full_jacobian); } + //(void)actor_state; // Calculate multiple scattering term // @TODO: Take stepping::config instead const auto joint_cov = stepping.calculate_ms_covariance( @@ -135,8 +136,8 @@ struct parameter_transporter : actor { // Add the multiple scattering term // new_cov = new_cov + stepping._joint_cov; new_cov = new_cov + joint_cov; - - /* + + /* std::cout << "Joint Cov" << std::endl; for (int i = 0; i < 6; i++) { for (int j = 0; j < 6; j++) { @@ -146,7 +147,7 @@ struct parameter_transporter : actor { } std::cout << std::endl; */ - + // Calculate surface-to-surface covariance transport stepping._bound_params.set_covariance(new_cov); } diff --git a/tests/integration_tests/cpu/propagator/jacobian_validation.cpp b/tests/integration_tests/cpu/propagator/jacobian_validation.cpp index 4b4819c8b..06dd8c267 100644 --- a/tests/integration_tests/cpu/propagator/jacobian_validation.cpp +++ b/tests/integration_tests/cpu/propagator/jacobian_validation.cpp @@ -341,7 +341,7 @@ struct bound_getter : actor { /// Numerically integrate the jacobian template bound_getter::state evaluate_bound_param( - const scalar detector_length, + const std::size_t mc_seed, const scalar detector_length, const bound_track_parameters& initial_param, const typename propagator_t::detector_type& det, const field_t& field, const scalar overstep_tolerance, const scalar on_surface_tolerance, @@ -371,6 +371,7 @@ bound_getter::state evaluate_bound_param( // Init propagator states for the reference track typename propagator_t::state state(initial_param, field, det); + state._stepping.rand_device.set_seed(mc_seed); // Run the propagation for the reference track state.do_debug = do_inspect; @@ -618,8 +619,8 @@ bound_track_parameters get_initial_parameter( template void evaluate_jacobian_difference( - const unsigned int trk_count, typename propagator_t::detector_type& det, - const scalar detector_length, + const std::size_t mc_seed, const unsigned int trk_count, + typename propagator_t::detector_type& det, const scalar detector_length, const bound_track_parameters& track, const field_t& field, const material volume_mat, const scalar overstep_tolerance, const scalar on_surface_tolerance, const scalar rk_tolerance, @@ -634,15 +635,15 @@ void evaluate_jacobian_difference( det.volumes()[0u].set_material(volume_mat); - const bool do_scatter = true; + const bool do_scatter = false; const bool do_covariance_transport = true; // No scattering for jacobian comparsion study // But perform covariance transport auto bound_getter = evaluate_bound_param( - detector_length, track, det, field, overstep_tolerance, + mc_seed, detector_length, track, det, field, overstep_tolerance, on_surface_tolerance, rk_tolerance, constraint_step, use_field_gradient, - !do_scatter, do_covariance_transport, do_inspect); + do_scatter, do_covariance_transport, do_inspect); const auto reference_param = bound_getter.m_param_departure; const auto final_param = bound_getter.m_param_destination; @@ -749,6 +750,7 @@ void evaluate_jacobian_difference( template void evaluate_covariance_transport( + const bool include_multiple_scattering, const std::size_t mc_seed, const unsigned int trk_count, typename propagator_t::detector_type& det, const scalar detector_length, const bound_track_parameters& track, const field_t& field, @@ -768,14 +770,14 @@ void evaluate_covariance_transport( track_copy.set_covariance(ini_cov); - const bool do_scatter = true; - const bool do_covariance_transport = true; + bool do_scatter = false; + bool do_covariance_transport = true; // Reference track (Without Scattering but with covariance transport) auto bound_getter = evaluate_bound_param( - detector_length, track_copy, det, field, overstep_tolerance, + mc_seed, detector_length, track_copy, det, field, overstep_tolerance, on_surface_tolerance, rk_tolerance, constraint_step, use_field_gradient, - !do_scatter, do_covariance_transport); + do_scatter, do_covariance_transport); const auto reference_param = bound_getter.m_param_departure; const auto ini_vec = reference_param.vector(); @@ -811,11 +813,18 @@ void evaluate_covariance_transport( auto smeared_track = track_copy; smeared_track.set_vector(smeared_ini_vec); + if (include_multiple_scattering) { + do_scatter = true; + } else { + do_scatter = false; + } + do_covariance_transport = false; + // Smeared track (With scattering and no covariance transport) auto smeared_bound_getter = evaluate_bound_param( - detector_length, smeared_track, det, field, overstep_tolerance, + mc_seed, detector_length, smeared_track, det, field, overstep_tolerance, on_surface_tolerance, rk_tolerance, constraint_step, use_field_gradient, - do_scatter, !do_covariance_transport); + do_scatter, do_covariance_transport); // Get smeared final bound vector bound_vector_type smeared_fin_vec = @@ -1343,6 +1352,9 @@ int main(int argc, char** argv) { "Set log10(max_rk_tolerance_in_mm)"); desc.add_options()("mc-seed", po::value()->default_value(0u), "Monte-Carlo seed"); + desc.add_options()("include-multiple-scattering", + po::value()->default_value(false), + "Include multiple scattering"); desc.add_options()("verbose-level", po::value()->default_value(1), "Verbose level"); @@ -1379,6 +1391,8 @@ int main(int argc, char** argv) { const scalar log10_max_rk_tolerance = vm["log10-max-rk-tolerance-mm"].as() * unit::mm; const std::size_t mc_seed = vm["mc-seed"].as(); + const bool include_multiple_scattering = + vm["include-multiple-scattering"].as(); const int verbose_lvl = vm["verbose-level"].as(); std::vector log10_tols; @@ -1672,11 +1686,11 @@ int main(int argc, char** argv) { // Rectangle Inhomogeneous field with Material evaluate_jacobian_difference( - track_count, rect_det, detector_length, rect_bparam, - inhom_bfield, volume_mat, overstep_tol, on_surface_tol, - std::pow(10.f, log10_tols[i]), constraint_step_size, - h_sizes_rect, rect_files[i], ref_rel_diff, true, - do_inspect); + mc_seed, track_count, rect_det, detector_length, + rect_bparam, inhom_bfield, volume_mat, overstep_tol, + on_surface_tol, std::pow(10.f, log10_tols[i]), + constraint_step_size, h_sizes_rect, rect_files[i], + ref_rel_diff, true, do_inspect); dqopdqop_rel_diffs_rect[i].push_back(ref_rel_diff); } @@ -1691,8 +1705,8 @@ int main(int argc, char** argv) { // Rect Const field evaluate_jacobian_difference( - track_count, rect_det, detector_length, rect_bparam, - const_bfield, vacuum(), overstep_tol, + mc_seed, track_count, rect_det, detector_length, + rect_bparam, const_bfield, vacuum(), overstep_tol, on_surface_tol, rk_tol, constraint_step_size, h_sizes_rect, const_rect_file, ref_rel_diff); @@ -1707,24 +1721,25 @@ int main(int argc, char** argv) { // Rect Inhomogeneous field evaluate_jacobian_difference( - track_count, rect_det, detector_length, rect_bparam, - inhom_bfield, vacuum(), overstep_tol, + mc_seed, track_count, rect_det, detector_length, + rect_bparam, inhom_bfield, vacuum(), overstep_tol, on_surface_tol, rk_tol, constraint_step_size, h_sizes_rect, inhom_rect_file, ref_rel_diff); // Rectangle Inhomogeneous field with Material evaluate_jacobian_difference( - track_count, rect_det, detector_length, rect_bparam, - inhom_bfield, volume_mat, overstep_tol, on_surface_tol, - rk_tol, constraint_step_size, h_sizes_rect, + mc_seed, track_count, rect_det, detector_length, + rect_bparam, inhom_bfield, volume_mat, overstep_tol, + on_surface_tol, rk_tol, constraint_step_size, h_sizes_rect, inhom_rect_material_file, ref_rel_diff); // Rectangle Inhomogeneous field with Material (Covariance // transport) evaluate_covariance_transport( - track_count, rect_det, detector_length, rect_bparam, - inhom_bfield, volume_mat, overstep_tol, on_surface_tol, - rk_tol, constraint_step_size, rect_cov_transport_file); + include_multiple_scattering, mc_seed, track_count, rect_det, + detector_length, rect_bparam, inhom_bfield, volume_mat, + overstep_tol, on_surface_tol, rk_tol, constraint_step_size, + rect_cov_transport_file); } } @@ -1750,11 +1765,11 @@ int main(int argc, char** argv) { for (std::size_t i = 0u; i < log10_tols.size(); i++) { // Wire Inhomogeneous field with Material evaluate_jacobian_difference( - track_count, wire_det, detector_length, wire_bparam, - inhom_bfield, volume_mat, overstep_tol, on_surface_tol, - std::pow(10.f, log10_tols[i]), constraint_step_size, - h_sizes_wire, wire_files[i], ref_rel_diff, true, - do_inspect); + mc_seed, track_count, wire_det, detector_length, + wire_bparam, inhom_bfield, volume_mat, overstep_tol, + on_surface_tol, std::pow(10.f, log10_tols[i]), + constraint_step_size, h_sizes_wire, wire_files[i], + ref_rel_diff, true, do_inspect); dqopdqop_rel_diffs_wire[i].push_back(ref_rel_diff); } @@ -1769,8 +1784,8 @@ int main(int argc, char** argv) { // Wire Const field evaluate_jacobian_difference( - track_count, wire_det, detector_length, wire_bparam, - const_bfield, vacuum(), overstep_tol, + mc_seed, track_count, wire_det, detector_length, + wire_bparam, const_bfield, vacuum(), overstep_tol, on_surface_tol, rk_tol, constraint_step_size, h_sizes_wire, const_wire_file, ref_rel_diff); @@ -1785,23 +1800,24 @@ int main(int argc, char** argv) { // Wire Inhomogeneous field evaluate_jacobian_difference( - track_count, wire_det, detector_length, wire_bparam, - inhom_bfield, vacuum(), overstep_tol, + mc_seed, track_count, wire_det, detector_length, + wire_bparam, inhom_bfield, vacuum(), overstep_tol, on_surface_tol, rk_tol, constraint_step_size, h_sizes_wire, inhom_wire_file, ref_rel_diff); // Wire Inhomogeneous field with Material evaluate_jacobian_difference( - track_count, wire_det, detector_length, wire_bparam, - inhom_bfield, volume_mat, overstep_tol, on_surface_tol, - rk_tol, constraint_step_size, h_sizes_wire, + mc_seed, track_count, wire_det, detector_length, + wire_bparam, inhom_bfield, volume_mat, overstep_tol, + on_surface_tol, rk_tol, constraint_step_size, h_sizes_wire, inhom_wire_material_file, ref_rel_diff); // Wire Inhomogeneous field with Material (Covariance transport) evaluate_covariance_transport( - track_count, wire_det, detector_length, wire_bparam, - inhom_bfield, volume_mat, overstep_tol, on_surface_tol, - rk_tol, constraint_step_size, wire_cov_transport_file); + include_multiple_scattering, mc_seed, track_count, wire_det, + detector_length, wire_bparam, inhom_bfield, volume_mat, + overstep_tol, on_surface_tol, rk_tol, constraint_step_size, + wire_cov_transport_file); } } } diff --git a/tests/scripts/run_jacobian_validation.sh b/tests/scripts/run_jacobian_validation.sh index 6eb66d7b6..a1954943d 100644 --- a/tests/scripts/run_jacobian_validation.sh +++ b/tests/scripts/run_jacobian_validation.sh @@ -27,6 +27,9 @@ log10_on_surface_tol=-3 # Monte-Carlo seed mc_seed=0 +# include multiple scattering +include_multiple_scattering=false + # Skip the first phase skip_first_phase=false @@ -36,7 +39,7 @@ skip_second_phase=false # Verbose level verbose_level=1 -while getopts "hd:n:t:p:q:c:i:s:f:r:v:" arg; do +while getopts "hd:n:t:p:q:c:i:s:m:f:r:v:" arg; do case $arg in h) echo "" @@ -51,6 +54,7 @@ while getopts "hd:n:t:p:q:c:i:s:f:r:v:" arg; do echo "-c " echo "-i " echo "-s " + echo "-m " echo "-f " echo "-r " echo "-v " @@ -90,6 +94,10 @@ while getopts "hd:n:t:p:q:c:i:s:f:r:v:" arg; do mc_seed=$OPTARG echo "Monte-Carlo seed: ${mc_seed}" ;; + m) + include_multiple_scattering=$OPTARG + echo "Include multiple scattering in the second phase: ${include_multiple_scattering}" + ;; f) skip_first_phase=$OPTARG echo "Skip the first phase: ${skip_first_phase}" @@ -141,6 +149,7 @@ if [ "$skip_first_phase" = false ] ; then --log10-helix-tolerance=${log10_helix_tol} \ --log10-on-surface-tolerance=${log10_on_surface_tol} \ --mc-seed=${mc_seed} \ + --include-multiple-scattering=${include_multiple_scattering} \ --verbose-level=${verbose_level}" ${command_rk_tolerance} & done @@ -170,7 +179,8 @@ if [ "$skip_second_phase" = false ] ; then --log10-rk-tolerance=${log10_rk_tol} \ --log10-helix-tolerance=${log10_helix_tol} \ --log10-on-surface-tolerance=${log10_on_surface_tol} \ - --mc-seed=${mc_seed} + --mc-seed=${mc_seed} \ + --include-multiple-scattering=${include_multiple_scattering} \ --verbose-level=${verbose_level}" ${command_jacobi_validation} & done @@ -235,7 +245,7 @@ if [ "$skip_first_phase" = false ] && [ "$skip_second_phase" = false ]; then elif [ "$skip_first_phase" = true ] && [ "$skip_second_phase" = false ]; then # Run covariance_validation.C - root -q -l ../../../tests/validation/root/covariance_validation.C+O + root ../../../tests/validation/root/covariance_validation.C+O elif [ "$skip_first_phase" = false ] && [ "$skip_second_phase" = true ]; then From 767b604527048e44aeeedad760f2451f1f8aaffd Mon Sep 17 00:00:00 2001 From: beomki-yeo Date: Sun, 18 Feb 2024 01:54:18 +0100 Subject: [PATCH 13/13] Backup --- .../actors/parameter_transporter.hpp | 2 +- core/include/detray/propagator/rk_stepper.ipp | 27 ++++++++++++++++--- 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/core/include/detray/propagator/actors/parameter_transporter.hpp b/core/include/detray/propagator/actors/parameter_transporter.hpp index a2340eee2..fab159e64 100644 --- a/core/include/detray/propagator/actors/parameter_transporter.hpp +++ b/core/include/detray/propagator/actors/parameter_transporter.hpp @@ -136,7 +136,7 @@ struct parameter_transporter : actor { // Add the multiple scattering term // new_cov = new_cov + stepping._joint_cov; new_cov = new_cov + joint_cov; - + /* std::cout << "Joint Cov" << std::endl; for (int i = 0; i < 6; i++) { diff --git a/core/include/detray/propagator/rk_stepper.ipp b/core/include/detray/propagator/rk_stepper.ipp index 110526fa0..d48b80650 100644 --- a/core/include/detray/propagator/rk_stepper.ipp +++ b/core/include/detray/propagator/rk_stepper.ipp @@ -446,11 +446,11 @@ DETRAY_HOST_DEVICE auto detray::rk_stepper< magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t, inspector_t, array_t>::state::calculate_ms_covariance(bool do_scatter, - bool do_covariance_transport) - const -> bound_matrix { + bool do_covariance_transport) const + -> bound_matrix { // Joint covariance - bound_matrix joint_cov = matrix_operator().template zero<6,6>(); + bound_matrix joint_cov = matrix_operator().template zero<6, 6>(); // Return if there is no material if (this->_mat == vacuum()) { @@ -527,6 +527,27 @@ DETRAY_HOST_DEVICE auto detray::rk_stepper< // Set the joint covariance matrix_operator().template set_block<4, 4>(joint_cov, E, 0, 0); + // Add qop variance + scalar_type sigma_qop = + interaction_type().compute_energy_loss_landau_sigma_QOverP( + L, this->_mat, this->_pdg, this->_mass, this->_track.qop(), + this->_track.charge()); + + // Flip the sign if the step length is negative + scalar_type var_qop = sigma_qop * sigma_qop; + if (L < 0) { + var_qop = -var_qop; + } + + /* + scalar_type qop1 = this->_qop_i; + scalar_type qop2 = this->_track.qop(); + + std::cout << L << " " << sigma_qop << " " << -1/qop1 << " " << -1/qop2 << std::endl; + */ + + getter::element(joint_cov, e_bound_qoverp, e_bound_qoverp) += var_qop; + return joint_cov; }