From cba726247f45e135220d2caae0a700b9b911be78 Mon Sep 17 00:00:00 2001 From: Oliver Lantwin Date: Thu, 9 Jul 2026 14:12:04 +0200 Subject: [PATCH 1/2] feat: type the Pythia8 boundary with canonical quantities Pythia8's native units (GeV, mm, mm/c production time) become explicit: beam energies are ship::Energy, lifetime thresholds and tProd are mm/c quantities, and the mm/c -> ns conversion derives from the exact SI definition of c instead of a hand-typed 299.792458. Fixed-target and Pythia8 source configuration is read via get_quantity, so each parameter's unit is defined in exactly one place. The only numeric change is the tProd conversion (division by 299.792458 becomes multiplication by the c-derived factor): at most 1 ulp on ~20% of values. Validation histograms of a 200-event pythia8_only run are bitwise identical to an origin/main build. --- CMakeLists.txt | 9 ++++--- benchmarks/pythia8_benchmark.cpp | 14 +++++++---- src/fixed_target_source.cpp | 39 +++++++++++++++++------------ src/pythia8_source.cpp | 8 +++--- src/pythia_common.hpp | 42 +++++++++++++++++++++----------- 5 files changed, 72 insertions(+), 40 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ea4b542..9831688 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/benchmarks/pythia8_benchmark.cpp b/benchmarks/pythia8_benchmark.cpp index aaadc8f..1101a30 100644 --- a/benchmarks/pythia8_benchmark.cpp +++ b/benchmarks/pythia8_benchmark.cpp @@ -78,7 +78,7 @@ void print_stats(char const* label, Stats const& s) { template 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(); diff --git a/src/fixed_target_source.cpp b/src/fixed_target_source.cpp index 644844f..b3a5318 100644 --- a/src/fixed_target_source.cpp +++ b/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 pythia_pp_; std::unique_ptr 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("beam_energy", 400.0); + auto beam_energy = + aegir::get_quantity(config, "beam_energy", 400.0 * su::GeV); auto target_z = config.get("target_z", 74); auto target_a = config.get("target_a", 184); - auto target_z_start = config.get("target_z_start", 0.0); - auto target_z_end = config.get("target_z_end", 1164.0); - auto interaction_length = config.get("interaction_length", 191.9); - auto tau0_threshold = config.get("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( "fixed_target", xml_dir, beam_energy, target_z, target_a, target_z_start, diff --git a/src/pythia8_source.cpp b/src/pythia8_source.cpp index 0e7ca2a..05af01f 100644 --- a/src/pythia8_source.cpp +++ b/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(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("beam_energy", 400.0); + auto beam_energy = + aegir::get_quantity(config, "beam_energy", 400.0 * ship::units::GeV); auto process = config.get("process", std::string{"SoftQCD:inelastic"}); auto parallel = config.get("parallel", false); diff --git a/src/pythia_common.hpp b/src/pythia_common.hpp index 0bb3641..8a22328 100644 --- a/src/pythia_common.hpp +++ b/src/pythia_common.hpp @@ -14,6 +14,7 @@ #include +#include #include #include #include @@ -21,26 +22,35 @@ namespace aegir { +// Pythia8's native units: energies/momenta in GeV, positions in mm, +// production times in mm/c. +using PythiaTime = mp_units::quantity; + // 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 -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 -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,7 +71,7 @@ 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 @@ -69,10 +79,12 @@ void next_event(Pythia& pythia, std::string_view source_name, // generally are not. This makes motherId a valid index into the emitted // collection rather than a dangling reference into the discarded record. template -std::vector extract_particles(Pythia8::Event const& event, - double z_offset = 0.0) { +std::vector extract_particles( + Pythia8::Event const& event, ship::Length z_offset = ship::Length::zero()) { + namespace su = ship::units; std::vector 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 out_index(static_cast(event.size()), -1); @@ -85,11 +97,13 @@ std::vector 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); } From 9218e28865132693ba9256588520cbc04de474de Mon Sep 17 00:00:00 2001 From: Oliver Lantwin Date: Thu, 9 Jul 2026 14:16:07 +0200 Subject: [PATCH 2/2] feat: type the GENIE and particle-gun boundaries GENIE's SI vertex (m, s) converts to the canonical mm/ns through quantities instead of hand-typed 1e3/1e9 factors (exact powers of ten, value-identical). Particle-gun configuration is read via get_quantity (momenta in GeV/c, angle in rad, vertex in mm); sampling stays on raw numbers at the RNG call so the Philox stream is bit-identical. Validation: determinism check passes; a 100-event gun_only run is bitwise identical to an origin/main build (order-independent event comparison; the parallel writer emits events in arbitrary order). --- CMakeLists.txt | 4 ++-- src/genie_reader_source.cpp | 14 +++++++++---- src/particle_gun_source.cpp | 39 +++++++++++++++++++++++-------------- 3 files changed, 36 insertions(+), 21 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9831688..f284a7f 100644 --- a/CMakeLists.txt +++ b/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} ) @@ -157,7 +157,7 @@ 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) diff --git a/src/genie_reader_source.cpp b/src/genie_reader_source.cpp index d3f9da0..a93e56b 100644 --- a/src/genie_reader_source.cpp +++ b/src/genie_reader_source.cpp @@ -15,6 +15,7 @@ #include #include +#include #include #include #include @@ -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 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 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 particles; particles.reserve(static_cast(n)); diff --git a/src/particle_gun_source.cpp b/src/particle_gun_source.cpp index 7b9fd41..d1b52ad 100644 --- a/src/particle_gun_source.cpp +++ b/src/particle_gun_source.cpp @@ -8,6 +8,7 @@ // Each event generates a single particle with fixed or randomised kinematics. #include +#include #include #include #include @@ -15,25 +16,31 @@ #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 vertex) + ParticleGun(int pdg, ship::Momentum p_min, ship::Momentum p_max, + ship::Angle max_theta, ship::Vec3 vertex) : pdg_{pdg}, p_min_{p_min}, p_max_{p_max}, max_theta_{max_theta}, - vertex_{vertex} {} + vertex_{ship::raw(vertex)} {} std::vector generate(phlex::data_cell_index const& id) { auto event_number = static_cast(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,8 +69,9 @@ class ParticleGun : public phlex::source { private: int pdg_; - double p_min_, p_max_, max_theta_; - std::array vertex_; + ship::Momentum p_min_, p_max_; + ship::Angle max_theta_; + std::array vertex_; // canonical mm }; } // namespace @@ -71,14 +79,15 @@ class ParticleGun : public phlex::source { PHLEX_REGISTER_SOURCE(s, config) { using namespace phlex; - auto pdg = config.get("pdg", 13); // muon - auto p_min = config.get("p_min", 10.0); // GeV - auto p_max = config.get("p_max", 100.0); - auto max_theta = config.get("max_theta", 0.1); // rad - auto vx = config.get("vertex_x", 0.0); - auto vy = config.get("vertex_y", 0.0); - auto vz = config.get("vertex_z", -500.0); // mm, upstream of target + auto pdg = config.get("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("particle_gun", pdg, p_min, p_max, max_theta, - std::array{vx, vy, vz}); + ship::Vec3{vx, vy, vz}); }