feat: add canonical unit vocabulary via mp-units

CMakeLists.txt

@@ -6,6 +6,7 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
 find_package(ROOT REQUIRED COMPONENTS Core RIO)
+find_package(mp-units 2.5 CONFIG REQUIRED)
 
 add_library(SHiPDataModel SHARED)
 target_sources(
@@ -24,6 +25,22 @@ target_sources(
 target_link_libraries(SHiPDataModel PUBLIC ROOT::Core)
 set_target_properties(SHiPDataModel PROPERTIES LINKER_LANGUAGE CXX)
 
+# Canonical unit vocabulary (mp-units quantity types + views over the data
+# classes). Headers-only and deliberately separate from SHiPDataModel: the
+# quantity types must never enter the ROOT dictionary (ROOT cannot stream
+# them), and units-only consumers must not pay for a ROOT linkage.
+add_library(SHiPUnits INTERFACE)
+add_library(SHiP::SHiPUnits ALIAS SHiPUnits)
+target_sources(
+    SHiPUnits
+    INTERFACE
+        FILE_SET headers
+            TYPE HEADERS
+            BASE_DIRS include
+            FILES include/SHiP/Units.hpp include/SHiP/QuantityView.hpp
+)
+target_link_libraries(SHiPUnits INTERFACE mp-units::mp-units)
+
 root_generate_dictionary(
   G__SHiPDataModel
   SHiP/MCParticle.hpp
@@ -45,6 +62,7 @@ install(
     LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
     FILE_SET headers
 )
+install(TARGETS SHiPUnits EXPORT SHiPDataModelTargets FILE_SET headers)
 install(
     FILES
         ${CMAKE_CURRENT_BINARY_DIR}/libSHiPDataModel_rdict.pcm

cmake/SHiPDataModelConfig.cmake.in

@@ -2,6 +2,7 @@
 
 include(CMakeFindDependencyMacro)
 find_dependency(ROOT COMPONENTS Core RIO)
+find_dependency(mp-units 2.5 CONFIG)
 
 include("${CMAKE_CURRENT_LIST_DIR}/SHiPDataModelTargets.cmake")
 

include/SHiP/QuantityView.hpp

@@ -0,0 +1,131 @@
+#pragma once
+
+#include "SHiP/MCParticle.hpp"
+#include "SHiP/RecParticle.hpp"
+#include "SHiP/SimHit.hpp"
+#include "SHiP/SimParticle.hpp"
+#include "SHiP/Units.hpp"
+
+/// Quantity views over the persisted data classes.
+///
+/// The data classes stay plain-double PODs (ROOT streaming cannot handle
+/// mp-units members); these free functions convert to and from the canonical
+/// quantity types so computation code never touches an implicit unit. All
+/// conversions are bitwise (the storage unit is the quantity's unit).
+namespace ship::view {
+
+// --- MCParticle -----------------------------------------------------------
+[[nodiscard]] inline Vec3<Length> vertex(SHiP::MCParticle const& p) {
+  return vecOf<Length>(p.vertex);
+}
+[[nodiscard]] inline Vec3<Momentum> momentum(SHiP::MCParticle const& p) {
+  return vecOf<Momentum>(p.momentum);
+}
+/// Total energy.
+[[nodiscard]] inline Energy energy(SHiP::MCParticle const& p) {
+  return quantityOf<Energy>(p.energy);
+}
+[[nodiscard]] inline Time time(SHiP::MCParticle const& p) {
+  return quantityOf<Time>(p.time);
+}
+inline void setVertex(SHiP::MCParticle& p, Vec3<Length> const& v) {
+  p.vertex = raw(v);
+}
+inline void setMomentum(SHiP::MCParticle& p, Vec3<Momentum> const& v) {
+  p.momentum = raw(v);
+}
+inline void setEnergy(SHiP::MCParticle& p, Energy e) { p.energy = raw(e); }
+inline void setTime(SHiP::MCParticle& p, Time t) { p.time = raw(t); }
+
+// --- SimHit ---------------------------------------------------------------
+[[nodiscard]] inline Vec3<Length> position(SHiP::SimHit const& h) {
+  return vecOf<Length>(h.position);
+}
+[[nodiscard]] inline Vec3<Momentum> momentum(SHiP::SimHit const& h) {
+  return vecOf<Momentum>(h.momentum);
+}
+[[nodiscard]] inline Energy energyDeposit(SHiP::SimHit const& h) {
+  return quantityOf<Energy>(h.energyDeposit);
+}
+[[nodiscard]] inline Time time(SHiP::SimHit const& h) {
+  return quantityOf<Time>(h.time);
+}
+[[nodiscard]] inline Length pathLength(SHiP::SimHit const& h) {
+  return quantityOf<Length>(h.pathLength);
+}
+inline void setPosition(SHiP::SimHit& h, Vec3<Length> const& v) {
+  h.position = raw(v);
+}
+inline void setMomentum(SHiP::SimHit& h, Vec3<Momentum> const& v) {
+  h.momentum = raw(v);
+}
+inline void setEnergyDeposit(SHiP::SimHit& h, Energy e) {
+  h.energyDeposit = raw(e);
+}
+inline void setTime(SHiP::SimHit& h, Time t) { h.time = raw(t); }
+inline void setPathLength(SHiP::SimHit& h, Length l) { h.pathLength = raw(l); }
+
+// --- SimParticle ----------------------------------------------------------
+[[nodiscard]] inline Vec3<Length> vertex(SHiP::SimParticle const& p) {
+  return vecOf<Length>(p.vertex);
+}
+[[nodiscard]] inline Vec3<Length> endpoint(SHiP::SimParticle const& p) {
+  return vecOf<Length>(p.endpoint);
+}
+[[nodiscard]] inline Vec3<Momentum> momentum(SHiP::SimParticle const& p) {
+  return vecOf<Momentum>(p.momentum);
+}
+/// Initial kinetic energy.
+[[nodiscard]] inline Energy energy(SHiP::SimParticle const& p) {
+  return quantityOf<Energy>(p.energy);
+}
+[[nodiscard]] inline Time time(SHiP::SimParticle const& p) {
+  return quantityOf<Time>(p.time);
+}
+inline void setVertex(SHiP::SimParticle& p, Vec3<Length> const& v) {
+  p.vertex = raw(v);
+}
+inline void setEndpoint(SHiP::SimParticle& p, Vec3<Length> const& v) {
+  p.endpoint = raw(v);
+}
+inline void setMomentum(SHiP::SimParticle& p, Vec3<Momentum> const& v) {
+  p.momentum = raw(v);
+}
+inline void setEnergy(SHiP::SimParticle& p, Energy e) { p.energy = raw(e); }
+inline void setTime(SHiP::SimParticle& p, Time t) { p.time = raw(t); }
+
+// --- RecParticle ----------------------------------------------------------
+[[nodiscard]] inline Vec3<Length> vertex(SHiP::RecParticle const& p) {
+  return vecOf<Length>(p.vertex);
+}
+[[nodiscard]] inline Vec3<Length> endpoint(SHiP::RecParticle const& p) {
+  return vecOf<Length>(p.endpoint);
+}
+[[nodiscard]] inline Vec3<Momentum> momentum(SHiP::RecParticle const& p) {
+  return vecOf<Momentum>(p.momentum);
+}
+/// Initial kinetic energy.
+[[nodiscard]] inline Energy energy(SHiP::RecParticle const& p) {
+  return quantityOf<Energy>(p.energy);
+}
+[[nodiscard]] inline Time time(SHiP::RecParticle const& p) {
+  return quantityOf<Time>(p.time);
+}
+/// Impact parameter wrt the primary vertex.
+[[nodiscard]] inline Length ipPV(SHiP::RecParticle const& p) {
+  return quantityOf<Length>(p.ipPV);
+}
+inline void setVertex(SHiP::RecParticle& p, Vec3<Length> const& v) {
+  p.vertex = raw(v);
+}
+inline void setEndpoint(SHiP::RecParticle& p, Vec3<Length> const& v) {
+  p.endpoint = raw(v);
+}
+inline void setMomentum(SHiP::RecParticle& p, Vec3<Momentum> const& v) {
+  p.momentum = raw(v);
+}
+inline void setEnergy(SHiP::RecParticle& p, Energy e) { p.energy = raw(e); }
+inline void setTime(SHiP::RecParticle& p, Time t) { p.time = raw(t); }
+inline void setIpPV(SHiP::RecParticle& p, Length l) { p.ipPV = raw(l); }
+
+}  // namespace ship::view

include/SHiP/Units.hpp

@@ -0,0 +1,92 @@
+#pragma once
+
+#include <mp-units/systems/hep.h>  // exports framework.h; includes systems/si.h
+
+#include <array>
+
+/// Canonical ShipSoft unit vocabulary.
+///
+/// Storage convention (matches the data-class doxygen contracts and Key4hep):
+/// positions and lengths [mm], momenta [GeV/c], energies [GeV], times [ns].
+/// Persisted classes keep plain doubles in these units; computation code
+/// works with the quantity types below and converts at the boundaries (see
+/// SHiP/QuantityView.hpp).
+///
+/// NEVER `using namespace ship::units` in a translation unit that includes
+/// CLHEP or Geant4 headers: G4SystemOfUnits.hh injects global `mm`, `GeV`,
+/// `tesla`, ... that collide with these names. Always qualify.
+namespace ship {
+
+namespace units {
+
+// Length
+inline constexpr auto mm = mp_units::si::milli<mp_units::si::metre>;
+inline constexpr auto cm = mp_units::si::centi<mp_units::si::metre>;
+inline constexpr auto m = mp_units::si::metre;
+
+// Time
+inline constexpr auto ns = mp_units::si::nano<mp_units::si::second>;
+inline constexpr auto s = mp_units::si::second;
+
+// Angle
+inline constexpr auto rad = mp_units::si::radian;
+
+// Energy
+inline constexpr auto MeV = mp_units::hep::unit_symbols::MeV;
+inline constexpr auto GeV = mp_units::hep::unit_symbols::GeV;
+
+// Speed of light as a named unit ("c", exact by SI definition). Building
+// momentum, mass and time-of-flight units from it keeps every conversion
+// factor derived rather than hand-typed.
+inline constexpr auto c = mp_units::hep::speed_of_light;
+
+// Momentum and mass
+inline constexpr auto GeV_per_c = GeV / c;
+inline constexpr auto GeV_per_c2 = GeV / mp_units::square(c);
+
+// Time of flight per unit length (Pythia8 production time, lifetime cuts)
+inline constexpr auto mm_per_c = mm / c;
+
+// Magnetic flux density
+inline constexpr auto tesla = mp_units::si::tesla;
+
+}  // namespace units
+
+/// Canonical quantity types: double representation, storage unit baked in.
+using Length = mp_units::quantity<units::mm, double>;
+using Time = mp_units::quantity<units::ns, double>;
+using Energy = mp_units::quantity<units::GeV, double>;
+using Momentum = mp_units::quantity<units::GeV_per_c, double>;
+using Mass = mp_units::quantity<units::GeV_per_c2, double>;
+using Angle = mp_units::quantity<units::rad, double>;
+using MagneticField = mp_units::quantity<units::tesla, double>;
+
+template <typename Q>
+using Vec3 = std::array<Q, 3>;
+
+/// Wrap a raw scalar, interpreting it in Q's canonical unit (bitwise).
+template <typename Q>
+[[nodiscard]] constexpr Q quantityOf(double v) {
+  return v * Q::reference;
+}
+
+/// Unwrap a scalar quantity to a raw double in Q's canonical unit (bitwise).
+template <typename Q>
+[[nodiscard]] constexpr double raw(Q q) {
+  return q.numerical_value_in(Q::unit);
+}
+
+/// Wrap a raw 3-vector, interpreting it in Q's canonical unit (bitwise).
+template <typename Q>
+[[nodiscard]] constexpr Vec3<Q> vecOf(std::array<double, 3> const& v) {
+  return {v[0] * Q::reference, v[1] * Q::reference, v[2] * Q::reference};
+}
+
+/// Unwrap a quantity 3-vector to raw doubles in Q's canonical unit (bitwise).
+template <typename Q>
+[[nodiscard]] constexpr std::array<double, 3> raw(Vec3<Q> const& v) {
+  return {v[0].numerical_value_in(Q::unit), v[1].numerical_value_in(Q::unit),
+          v[2].numerical_value_in(Q::unit)};
+}
+
+}  // namespace ship

pixi.toml

@@ -16,6 +16,8 @@ ninja = "*"
 git = "*"
 root = ">=6.40"
 root_cxx_standard = "==23"
+# Required by the SHiP::SHiPUnits target (find_package(mp-units CONFIG)).
+mp-units = ">=2.5"
 
 [tasks]
 configure = "cmake -S . -B build -G Ninja -DCMAKE_BUILD_TYPE=RelWithDebInfo"

tests/CMakeLists.txt

@@ -18,3 +18,8 @@ target_link_libraries(
     PRIVATE ${link_dictionary} ROOT::Core ROOT::RIO ROOT::ROOTNTuple
 )
 add_test(NAME rntuple_io COMMAND test_rntuple_io)
+
+# Units and quantity views are headers-only: no dictionary, no ROOT.
+add_executable(test_units test_units.cpp)
+target_link_libraries(test_units PRIVATE SHiP::SHiPUnits)
+add_test(NAME units COMMAND test_units)