feat: type the GENIE and particle-gun boundaries

CMakeLists.txt

@@ -85,7 +85,7 @@ endfunction()
 
 aegir_add_plugin(
     particle_gun_source
-    LIBS SHiP::SHiPDataModel
+    LIBS SHiP::SHiPDataModel SHiP::SHiPUnits
     INCLUDES ${Random123_INCLUDE_DIR}
 )
 
@@ -140,12 +140,12 @@ aegir_add_plugin(field_null_provider LIBS SHiPFieldService::SHiPFieldService)
 
 aegir_add_plugin(
     pythia8_source
-    LIBS SHiP::SHiPDataModel Pythia8::Pythia8
+    LIBS SHiP::SHiPDataModel SHiP::SHiPUnits Pythia8::Pythia8
 )
 
 aegir_add_plugin(
     fixed_target_source
-    LIBS SHiP::SHiPDataModel Pythia8::Pythia8
+    LIBS SHiP::SHiPDataModel SHiP::SHiPUnits Pythia8::Pythia8
     INCLUDES ${Random123_INCLUDE_DIR}
 )
 
@@ -157,12 +157,15 @@ aegir_add_plugin(
 
 aegir_add_plugin(
     genie_reader_source
-    LIBS SHiP::SHiPDataModel ROOT::RIO ROOT::Tree
+    LIBS SHiP::SHiPDataModel SHiP::SHiPUnits ROOT::RIO ROOT::Tree
 )
 
 # Standalone Pythia8 benchmark (no Phlex dependency)
 add_executable(pythia8_benchmark benchmarks/pythia8_benchmark.cpp)
-target_link_libraries(pythia8_benchmark PRIVATE Pythia8::Pythia8)
+target_link_libraries(
+    pythia8_benchmark
+    PRIVATE Pythia8::Pythia8 SHiP::SHiPUnits
+)
 target_include_directories(
     pythia8_benchmark
     PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/src

benchmarks/pythia8_benchmark.cpp

@@ -78,7 +78,7 @@ void print_stats(char const* label, Stats const& s) {
 template <typename T>
 void configure_pythia(T& pythia, double beam_energy, std::string const& process,
                       bool fairship) {
-  aegir::configure_beams(pythia, 2212, 2212, beam_energy);
+  aegir::configure_beams(pythia, 2212, 2212, beam_energy * ship::units::GeV);
   pythia.readString("Print:quiet = on");
 
   if (fairship) {
@@ -166,7 +166,8 @@ int main(int argc, char* argv[]) {
     Pythia8::Pythia pythia(xml_dir(), false);
     configure_pythia(pythia, cfg.beam_energy, cfg.process, cfg.fairship);
     pythia.init();
-    if (cfg.fairship) aegir::stabilise_long_lived(pythia, 1.0);
+    if (cfg.fairship)
+      aegir::stabilise_long_lived(pythia, 1.0 * ship::units::mm_per_c);
 
     // Warmup
     for (int i = 0; i < cfg.warmup; ++i) pythia.next();
@@ -187,7 +188,8 @@ int main(int argc, char* argv[]) {
     Pythia8::Pythia pythia(xml_dir(), false);
     configure_pythia(pythia, cfg.beam_energy, cfg.process, cfg.fairship);
     pythia.init();
-    if (cfg.fairship) aegir::stabilise_long_lived(pythia, 1.0);
+    if (cfg.fairship)
+      aegir::stabilise_long_lived(pythia, 1.0 * ship::units::mm_per_c);
 
     // Warmup
     for (int i = 0; i < cfg.warmup; ++i) {
@@ -216,7 +218,8 @@ int main(int argc, char* argv[]) {
     pythia.readString("Parallelism:numThreads = " +
                       std::to_string(cfg.threads));
     pythia.init();
-    if (cfg.fairship) aegir::stabilise_long_lived(pythia, 1.0);
+    if (cfg.fairship)
+      aegir::stabilise_long_lived(pythia, 1.0 * ship::units::mm_per_c);
 
     // Warmup
     int warmup_count = 0;
@@ -231,7 +234,8 @@ int main(int argc, char* argv[]) {
     pythia2.readString("Parallelism:numThreads = " +
                        std::to_string(cfg.threads));
     pythia2.init();
-    if (cfg.fairship) aegir::stabilise_long_lived(pythia2, 1.0);
+    if (cfg.fairship)
+      aegir::stabilise_long_lived(pythia2, 1.0 * ship::units::mm_per_c);
 
     int event_count = 0;
     auto t0 = Clock::now();

src/fixed_target_source.cpp

@@ -22,9 +22,12 @@
 #include "mc_particle_source.hpp"
 #include "philox_rng.hpp"
 #include "pythia_common.hpp"
+#include "units/config_units.hpp"
 
 namespace {
 
+namespace su = ship::units;
+
 void configure_processes(Pythia8::Pythia& pythia) {
   pythia.readString("SoftQCD:inelastic = on");
   pythia.readString("PhotonCollision:gmgm2mumu = on");
@@ -35,10 +38,10 @@ void configure_processes(Pythia8::Pythia& pythia) {
 
 class FixedTargetSource : public phlex::source {
  public:
-  FixedTargetSource(std::string const& xml_dir, double beam_energy,
-                    int target_z, int target_a, double target_z_start,
-                    double target_z_end, double interaction_length,
-                    double tau0_threshold)
+  FixedTargetSource(std::string const& xml_dir, ship::Energy beam_energy,
+                    int target_z, int target_a, ship::Length target_z_start,
+                    ship::Length target_z_end, ship::Length interaction_length,
+                    aegir::PythiaTime tau0_threshold)
       : target_z_{target_z},
         target_a_{target_a},
         target_z_start_{target_z_start},
@@ -73,10 +76,11 @@ class FixedTargetSource : public phlex::source {
     bool proton_target = rng.uniform() < z_over_a;
 
     // Sample interaction point z from truncated exponential
-    double target_length = target_z_end_ - target_z_start_;
+    ship::Length target_length = target_z_end_ - target_z_start_;
     double u = rng.uniform();
-    double exp_ratio = std::exp(-target_length / interaction_length_);
-    double z_interaction =
+    double exp_ratio = std::exp(-(target_length / interaction_length_)
+                                     .numerical_value_in(mp_units::one));
+    ship::Length z_interaction =
         target_z_start_ -
         interaction_length_ * std::log(1.0 - u * (1.0 - exp_ratio));
 
@@ -101,9 +105,9 @@ class FixedTargetSource : public phlex::source {
  private:
   int target_z_;
   int target_a_;
-  double target_z_start_;
-  double target_z_end_;
-  double interaction_length_;
+  ship::Length target_z_start_;
+  ship::Length target_z_end_;
+  ship::Length interaction_length_;
   std::unique_ptr<Pythia8::Pythia> pythia_pp_;
   std::unique_ptr<Pythia8::Pythia> pythia_pn_;
 };
@@ -117,13 +121,18 @@ PHLEX_REGISTER_SOURCE(s, config) {
     if (auto const* env = std::getenv("PYTHIA8DATA")) return std::string{env};
     return std::string{"../share/Pythia8/xmldoc"};
   }());
-  auto beam_energy = config.get<double>("beam_energy", 400.0);
+  auto beam_energy =
+      aegir::get_quantity(config, "beam_energy", 400.0 * su::GeV);
   auto target_z = config.get<int>("target_z", 74);
   auto target_a = config.get<int>("target_a", 184);
-  auto target_z_start = config.get<double>("target_z_start", 0.0);
-  auto target_z_end = config.get<double>("target_z_end", 1164.0);
-  auto interaction_length = config.get<double>("interaction_length", 191.9);
-  auto tau0_threshold = config.get<double>("tau0_threshold", 1.0);
+  auto target_z_start =
+      aegir::get_quantity(config, "target_z_start", 0.0 * su::mm);
+  auto target_z_end =
+      aegir::get_quantity(config, "target_z_end", 1164.0 * su::mm);
+  auto interaction_length =
+      aegir::get_quantity(config, "interaction_length", 191.9 * su::mm);
+  auto tau0_threshold =
+      aegir::get_quantity(config, "tau0_threshold", 1.0 * su::mm_per_c);
 
   s.add_source<FixedTargetSource>(
       "fixed_target", xml_dir, beam_energy, target_z, target_a, target_z_start,

src/genie_reader_source.cpp

@@ -15,6 +15,7 @@
 #include <TTree.h>
 
 #include <SHiP/MCParticle.hpp>
+#include <SHiP/Units.hpp>
 #include <algorithm>
 #include <memory>
 #include <stdexcept>
@@ -25,6 +26,8 @@
 
 namespace {
 
+namespace su = ship::units;
+
 // GENIE GHepStatus: 1 = stable final state (trackable by Geant4).
 constexpr int kStableFinalState = 1;
 
@@ -89,10 +92,13 @@ class GenieReaderSource : public phlex::source {
 
     // The interaction vertex is per event (EvtVtx, SI units); the per-particle
     // StdHepX4 positions are nuclear-scale offsets and irrelevant for
-    // tracking.
-    std::array<double, 3> const vertex{vtx_[0] * 1e3, vtx_[1] * 1e3,
-                                       vtx_[2] * 1e3};  // m -> mm
-    double const time = vtx_[3] * 1e9;                  // s -> ns
+    // tracking. GENIE's SI values acquire their unit here and convert to the
+    // canonical storage units (exact powers of ten).
+    std::array<double, 3> const vertex{
+        (vtx_[0] * su::m).numerical_value_in(su::mm),
+        (vtx_[1] * su::m).numerical_value_in(su::mm),
+        (vtx_[2] * su::m).numerical_value_in(su::mm)};
+    double const time = (vtx_[3] * su::s).numerical_value_in(su::ns);
 
     std::vector<SHiP::MCParticle> particles;
     particles.reserve(static_cast<std::size_t>(n));

src/particle_gun_source.cpp

@@ -8,32 +8,39 @@
 // Each event generates a single particle with fixed or randomised kinematics.
 
 #include <SHiP/MCParticle.hpp>
+#include <SHiP/QuantityView.hpp>
 #include <cmath>
 #include <cstdint>
 #include <numbers>
 #include <vector>
 
 #include "mc_particle_source.hpp"
 #include "philox_rng.hpp"
+#include "units/config_units.hpp"
 
 namespace {
 
+namespace su = ship::units;
+
 class ParticleGun : public phlex::source {
  public:
-  ParticleGun(int pdg, double p_min, double p_max, double max_theta,
-              std::array<double, 3> vertex)
+  ParticleGun(int pdg, ship::Momentum p_min, ship::Momentum p_max,
+              ship::Angle max_theta, ship::Vec3<ship::Length> vertex)
       : pdg_{pdg},
         p_min_{p_min},
         p_max_{p_max},
         max_theta_{max_theta},
-        vertex_{vertex} {}
+        vertex_{ship::raw(vertex)} {}
 
   std::vector<SHiP::MCParticle> generate(phlex::data_cell_index const& id) {
     auto event_number = static_cast<std::uint32_t>(id.number());
     aegir::PhiloxRng rng{event_number};
 
-    double p = rng.uniform(p_min_, p_max_);
-    double theta = rng.uniform(0.0, max_theta_);
+    // Sampling stays on raw numbers so the Philox stream is bit-identical;
+    // the bounds are unwrapped in the canonical units on the same lines.
+    double p = rng.uniform(p_min_.numerical_value_in(su::GeV_per_c),
+                           p_max_.numerical_value_in(su::GeV_per_c));
+    double theta = rng.uniform(0.0, max_theta_.numerical_value_in(su::rad));
     double phi = rng.uniform(0.0, 2.0 * std::numbers::pi);
 
     SHiP::MCParticle mc;
@@ -62,23 +69,25 @@ class ParticleGun : public phlex::source {
 
  private:
   int pdg_;
-  double p_min_, p_max_, max_theta_;
-  std::array<double, 3> vertex_;
+  ship::Momentum p_min_, p_max_;
+  ship::Angle max_theta_;
+  std::array<double, 3> vertex_;  // canonical mm
 };
 
 }  // namespace
 
 PHLEX_REGISTER_SOURCE(s, config) {
   using namespace phlex;
 
-  auto pdg = config.get<int>("pdg", 13);           // muon
-  auto p_min = config.get<double>("p_min", 10.0);  // GeV
-  auto p_max = config.get<double>("p_max", 100.0);
-  auto max_theta = config.get<double>("max_theta", 0.1);  // rad
-  auto vx = config.get<double>("vertex_x", 0.0);
-  auto vy = config.get<double>("vertex_y", 0.0);
-  auto vz = config.get<double>("vertex_z", -500.0);  // mm, upstream of target
+  auto pdg = config.get<int>("pdg", 13);  // muon
+  auto p_min = aegir::get_quantity(config, "p_min", 10.0 * su::GeV_per_c);
+  auto p_max = aegir::get_quantity(config, "p_max", 100.0 * su::GeV_per_c);
+  auto max_theta = aegir::get_quantity(config, "max_theta", 0.1 * su::rad);
+  auto vx = aegir::get_quantity(config, "vertex_x", 0.0 * su::mm);
+  auto vy = aegir::get_quantity(config, "vertex_y", 0.0 * su::mm);
+  // Default: upstream of the target.
+  auto vz = aegir::get_quantity(config, "vertex_z", -500.0 * su::mm);
 
   s.add_source<ParticleGun>("particle_gun", pdg, p_min, p_max, max_theta,
-                            std::array<double, 3>{vx, vy, vz});
+                            ship::Vec3<ship::Length>{vx, vy, vz});
 }

src/pythia8_source.cpp

@@ -26,6 +26,7 @@
 
 #include "mc_particle_source.hpp"
 #include "pythia_common.hpp"
+#include "units/config_units.hpp"
 
 namespace {
 
@@ -35,7 +36,7 @@ namespace {
 
 class Pythia8Source : public phlex::source {
  public:
-  Pythia8Source(std::string const& xml_dir, double beam_energy,
+  Pythia8Source(std::string const& xml_dir, ship::Energy beam_energy,
                 std::string const& process) {
     pythia_ = std::make_unique<Pythia8::Pythia>(xml_dir, false);
     aegir::configure_beams(*pythia_, 2212, 2212, beam_energy);
@@ -69,7 +70,7 @@ class Pythia8Source : public phlex::source {
 
 class Pythia8MTSource : public phlex::source {
  public:
-  Pythia8MTSource(std::string const& xml_dir, double beam_energy,
+  Pythia8MTSource(std::string const& xml_dir, ship::Energy beam_energy,
                   std::string const& process, int num_threads, long num_events,
                   std::size_t max_queue_size)
       : max_queue_size_{max_queue_size} {
@@ -202,7 +203,8 @@ PHLEX_REGISTER_SOURCE(s, config) {
     if (auto const* env = std::getenv("PYTHIA8DATA")) return std::string{env};
     return std::string{"../share/Pythia8/xmldoc"};
   }());
-  auto beam_energy = config.get<double>("beam_energy", 400.0);
+  auto beam_energy =
+      aegir::get_quantity(config, "beam_energy", 400.0 * ship::units::GeV);
   auto process =
       config.get<std::string>("process", std::string{"SoftQCD:inelastic"});
   auto parallel = config.get<bool>("parallel", false);