feat: type the Pythia8 boundary with canonical quantities

CMakeLists.txt

@@ -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}
 )
 
@@ -162,7 +162,10 @@ aegir_add_plugin(
 
 # 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/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);

src/pythia_common.hpp

@@ -14,33 +14,43 @@
 
 #include <Pythia8/Pythia.h>
 
+#include <SHiP/Units.hpp>
 #include <stdexcept>
 #include <string>
 #include <string_view>
 #include <vector>
 
 namespace aegir {
 
+// Pythia8's native units: energies/momenta in GeV, positions in mm,
+// production times in mm/c.
+using PythiaTime = mp_units::quantity<ship::units::mm_per_c, double>;
+
 // Configure a fixed-target beam: beam A on a stationary target B (frameType 2,
 // eB = 0). Templated so it works for both Pythia8::Pythia and PythiaParallel.
 template <typename Pythia>
-void configure_beams(Pythia& pythia, int idA, int idB, double beam_energy) {
+void configure_beams(Pythia& pythia, int idA, int idB,
+                     ship::Energy beam_energy) {
   pythia.readString("Beams:idA = " + std::to_string(idA));
   pythia.readString("Beams:idB = " + std::to_string(idB));
   pythia.readString("Beams:frameType = 2");
-  pythia.readString("Beams:eA = " + std::to_string(beam_energy));
+  pythia.readString(
+      "Beams:eA = " +
+      std::to_string(beam_energy.numerical_value_in(ship::units::GeV)));
   pythia.readString("Beams:eB = 0.");
 }
 
-// Make long-lived particles (tau0 above threshold, in mm/c) stable so a
-// downstream simulation (e.g. Geant4) handles their decay. Guards against null
+// Make long-lived particles (tau0 above threshold) stable so a downstream
+// simulation (e.g. Geant4) handles their decay. Guards against null
 // particleData entries.
 template <typename Pythia>
-void stabilise_long_lived(Pythia& pythia, double tau0_threshold) {
+void stabilise_long_lived(Pythia& pythia, PythiaTime tau0_threshold) {
+  double const threshold =
+      tau0_threshold.numerical_value_in(ship::units::mm_per_c);
   for (auto it = pythia.particleData.begin(); it != pythia.particleData.end();
        ++it) {
     auto& entry = it->second;  // ParticleDataEntryPtr (shared_ptr-like)
-    if (entry && entry->tau0() > tau0_threshold) entry->setMayDecay(false);
+    if (entry && entry->tau0() > threshold) entry->setMayDecay(false);
   }
 }
 
@@ -61,18 +71,20 @@ void next_event(Pythia& pythia, std::string_view source_name,
 
 // Extract final-state particles from a Pythia event record into a vector of
 // MCParticle (any type exposing pdgCode/vertex/momentum/energy/time/motherId/
-// status). vertex z is shifted by z_offset (mm).
+// status). vertex z is shifted by z_offset.
 //
 // motherId is remapped from the full Pythia-record index to the index within
 // the returned vector, or -1 when the mother was not itself written out — the
 // common case, since only final-state particles are kept and their mothers
 // generally are not. This makes motherId a valid index into the emitted
 // collection rather than a dangling reference into the discarded record.
 template <typename MCParticle>
-std::vector<MCParticle> extract_particles(Pythia8::Event const& event,
-                                          double z_offset = 0.0) {
+std::vector<MCParticle> extract_particles(
+    Pythia8::Event const& event, ship::Length z_offset = ship::Length::zero()) {
+  namespace su = ship::units;
   std::vector<MCParticle> particles;
   particles.reserve(event.size());
+  double const z_offset_mm = z_offset.numerical_value_in(su::mm);
 
   // Pythia-record index -> output index for written (final-state) particles.
   std::vector<int> out_index(static_cast<std::size_t>(event.size()), -1);
@@ -85,11 +97,13 @@ std::vector<MCParticle> extract_particles(Pythia8::Event const& event,
 
     MCParticle mc;
     mc.pdgCode = p.id();
-    mc.vertex = {p.xProd(), p.yProd(), p.zProd() + z_offset};  // mm
-    mc.momentum = {p.px(), p.py(), p.pz()};                    // GeV
-    mc.energy = p.e();
-    mc.time = p.tProd() / 299.792458;  // mm/c -> ns
-    mc.motherId = p.mother1();         // record index, remapped below
+    // Pythia positions and momenta are already in the canonical units.
+    mc.vertex = {p.xProd(), p.yProd(), p.zProd() + z_offset_mm};  // mm
+    mc.momentum = {p.px(), p.py(), p.pz()};                       // GeV/c
+    mc.energy = p.e();                                            // GeV
+    // mm/c -> ns via the exact definition of c (no hand-typed constant).
+    mc.time = (p.tProd() * su::mm_per_c).numerical_value_in(su::ns);
+    mc.motherId = p.mother1();  // record index, remapped below
     mc.status = p.statusHepMC();
     particles.push_back(mc);
   }