traccc/simulation/simulator.hpp

0001 /** TRACCC library, part of the ACTS project (R&D line)
0002  *
0003  * (c) 2022-2026 CERN for the benefit of the ACTS project
0004  *
0005  * Mozilla Public License Version 2.0
0006  */
0007
0008 #pragma once
0009
0010 // Project include(s).
0011 #include "traccc/simulation/smearing_writer.hpp"
0012 #include "traccc/utils/logging.hpp"
0013 #include "traccc/utils/particle.hpp"
0014 #include "traccc/utils/propagation.hpp"
0015
0016 // Detray include(s).
0017 #include <detray/test/utils/random_scatterer.hpp>
0018
0019 // System include(s).
0020 #include <limits>
0021 #include <memory>
0022
0023 namespace traccc {
0024
0025 template <typename detector_t, typename bfield_t, typename track_generator_t,
0026           typename writer_t>
0027 struct simulator {
0028   using algebra_type = typename detector_t::algebra_type;
0029   using scalar_type = typename detector_t::scalar_type;
0030   using bfield_type = bfield_t;
0031
0032   struct config {
0033     detray::propagation::config propagation;
0034
0035     /// Particle hypothesis
0036     traccc::pdg_particle<scalar_type> ptc_type{traccc::muon<scalar_type>()};
0037
0038     // Simulation setup
0039     bool do_energy_loss = true;
0040     bool do_multiple_scattering = true;
0041     bool m_is_min_pT = false;
0042     scalar_type m_min_p = 10.f * traccc::unit<scalar_type>::MeV;
0043
0044     /// Set the momentum limit to @param p
0045     inline void min_p(const scalar_type p) {
0046       m_is_min_pT = false;
0047       m_min_p = p;
0048     }
0049
0050     /// Set the transverse momentum limit to @param p
0051     inline void min_pT(const scalar_type p) {
0052       m_is_min_pT = true;
0053       m_min_p = p;
0054     }
0055   };
0056
0057   using actor_chain_type = detray::actor_chain<
0058       detray::actor::parameter_updater<
0059           algebra_type, detray::actor::random_scatterer<algebra_type>,
0060           writer_t>,
0061       detray::actor::momentum_aborter<scalar_type>>;
0062
0063   using navigator_type = detray::caching_navigator<detector_t>;
0064   using stepper_type = detray::rk_stepper<
0065       typename bfield_type::view_t, algebra_type,
0066       detray::constrained_step<detray::dscalar<algebra_type>>>;
0067   using propagator_type =
0068       detray::propagator<stepper_type, navigator_type, actor_chain_type>;
0069
0070   simulator(const detray::pdg_particle<scalar>& ptc, std::size_t events,
0071             const detector_t& det, const bfield_type& field,
0072             track_generator_t&& track_gen,
0073             typename writer_t::config&& writer_cfg,
0074             const std::string directory = "")
0075       : m_events(events),
0076         m_directory(directory),
0077         m_detector(det),
0078         m_field(field),
0079         m_track_generator(
0080             std::make_unique<track_generator_t>(std::move(track_gen))),
0081         m_writer_cfg(writer_cfg) {
0082     m_cfg.ptc_type = ptc;
0083     m_track_generator->config().charge(ptc.charge());
0084
0085     // Turn off tracking features
0086     m_cfg.propagation.stepping.do_covariance_transport = false;
0087     m_cfg.propagation.stepping.use_eloss_gradient = false;
0088     m_cfg.propagation.stepping.use_field_gradient = false;
0089     m_cfg.propagation.navigation.estimate_scattering_noise = false;
0090
0091     m_updater_state =
0092         detray::actor::parameter_updater_state<algebra_type>{m_cfg.propagation};
0093   }
0094
0095   config& get_config() { return m_cfg; }
0096
0097   void run() {
0098     TRACCC_VERBOSE_HOST("Running fast simulation...");
0099
0100     // Update the actor config
0101     if (m_cfg.m_is_min_pT) {
0102       m_aborter_state.min_pT(m_cfg.m_min_p);
0103       m_aborter_state.min_p(0.f);
0104     } else {
0105       m_aborter_state.min_p(m_cfg.m_min_p);
0106       m_aborter_state.min_pT(0.f);
0107     }
0108     m_scatterer.do_energy_loss = m_cfg.do_energy_loss;
0109     m_scatterer.do_multiple_scattering = m_cfg.do_multiple_scattering;
0110
0111     for (std::size_t event_id = 0u; event_id < m_events; event_id++) {
0112       typename writer_t::state writer_state(event_id, std::move(m_writer_cfg),
0113                                             m_directory);
0114
0115       // Set random seed
0116       m_scatterer.set_seed(event_id);
0117       writer_state.set_seed(event_id);
0118
0119       auto actor_states = detray::tie(m_updater_state, m_scatterer,
0120                                       m_aborter_state, writer_state);
0121
0122       for (auto track : *m_track_generator.get()) {
0123         m_updater_state.init(track);
0124
0125         writer_state.write_particle(
0126             track, detail::correct_particle_hypothesis(m_cfg.ptc_type, track));
0127
0128         typename propagator_type::state propagation(track, m_field, m_detector);
0129         propagation.set_particle(
0130             detail::correct_particle_hypothesis(m_cfg.ptc_type, track));
0131
0132         propagator_type p(m_cfg.propagation);
0133
0134         // Set overstep tolerance and stepper constraint
0135         propagation.stepping()
0136             .template set_constraint<detray::step::constraint::e_accuracy>(
0137                 m_cfg.propagation.stepping.step_constraint);
0138
0139         p.propagate(propagation, actor_states);
0140
0141         // Increase the particle id
0142         writer_state.particle_id++;
0143       }
0144     }
0145   }
0146
0147  private:
0148   config m_cfg;
0149   std::size_t m_events{0u};
0150   std::string m_directory = "";
0151   const detector_t& m_detector;
0152   const typename bfield_type::view_t m_field;
0153   std::unique_ptr<track_generator_t> m_track_generator;
0154   typename writer_t::config m_writer_cfg;
0155
0156   /// Actor states
0157   typename detray::actor::momentum_aborter<scalar_type>::state
0158       m_aborter_state{};
0159   typename detray::actor::parameter_updater<algebra_type>::state
0160       m_updater_state{};
0161   typename detray::actor::random_scatterer<algebra_type>::state m_scatterer{};
0162 };
0163
0164 }  // namespace traccc