From 1e035cda2c73311f3854c3b78173ee1aff9d7704 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Tue, 5 May 2026 10:45:11 +0200 Subject: [PATCH 01/28] Include new C++ classes in opengate_lib for field management. This: - Allows fields to be defined relative to the local coordinate system of the physical volumes. - Paves the way for the mapped field implementation --- core/opengate_core/opengate_core.cpp | 6 + core/opengate_core/opengate_lib/GateField.cpp | 132 ++++++++++++++++++ core/opengate_core/opengate_lib/GateField.h | 69 +++++++++ .../opengate_lib/GateMagneticField.cpp | 50 +++++++ .../opengate_lib/GateMagneticField.h | 58 ++++++++ .../opengate_lib/pyGateMagneticField.cpp | 39 ++++++ 6 files changed, 354 insertions(+) create mode 100644 core/opengate_core/opengate_lib/GateField.cpp create mode 100644 core/opengate_core/opengate_lib/GateField.h create mode 100644 core/opengate_core/opengate_lib/GateMagneticField.cpp create mode 100644 core/opengate_core/opengate_lib/GateMagneticField.h create mode 100644 core/opengate_core/opengate_lib/pyGateMagneticField.cpp diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index 8eb2cfa39b..f09ff7a458 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -210,6 +210,10 @@ void init_G4MagInt_Driver(py::module &); void init_G4ChordFinder(py::module &); +void init_GateMagneticField(py::module &); + +// void init_GateMappedMagneticField(py::module &); + // geometry/solids void init_G4Box(py::module &); @@ -584,6 +588,8 @@ PYBIND11_MODULE(opengate_core, m) { init_G4VIntegrationDriver(m); init_G4MagInt_Driver(m); init_G4ChordFinder(m); + init_GateMagneticField(m); + // init_GateMappedMagneticField(m); init_G4Box(m); init_G4Ellipsoid(m); diff --git a/core/opengate_core/opengate_lib/GateField.cpp b/core/opengate_core/opengate_lib/GateField.cpp new file mode 100644 index 0000000000..eb8367b8d9 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateField.cpp @@ -0,0 +1,132 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include "GateField.h" + +#include "G4VSolid.hh" +#include "globals.hh" // G4Exception + +#include +#include +#include +#include +#include +#include + +// constructor +GateField::GateField( + const G4VSolid *solid, + std::vector translations, + std::vector rotations, + double deltaChordMM +) : m_solid(solid), + m_fallbackFatalDistanceMM(5.0 * std::sqrt(8.0 * kMaxCurvatureRadiusMM * deltaChordMM)) +{ + // sanity-check the inputs before caching + if (solid == nullptr) + throw std::invalid_argument("GateField: solid must not be null"); + + if (translations.size() != rotations.size() || translations.empty()) + throw std::invalid_argument("GateField: translations and rotations must be " + "non-empty and have the same size"); + + + // initial cache the world-to-local transforms for every physical placement of the physical volume + m_transforms.reserve(translations.size()); + for (std::size_t i = 0; i < translations.size(); ++i) + m_transforms.emplace_back(rotations[i].inverse(), translations[i]); + +} + +// recache the world-to-local transforms (e.g. after a geometry change between runs) +void GateField::SetTransforms( + std::vector translations, + std::vector rotations +) { + + // sanity-check the inputs before caching + if (translations.size() != rotations.size() || translations.empty()) + throw std::invalid_argument( + "GateField::SetTransforms: translations and rotations must be " + "non-empty and have the same size"); + + // recache the world-to-local transforms for every physical placement of the logical volume + m_transforms.clear(); + m_transforms.reserve(translations.size()); + for (std::size_t i = 0; i < translations.size(); ++i) + m_transforms.emplace_back(rotations[i].inverse(), translations[i]); + +} + +// find the local coordinates of worldPoint in the containing placement of the +// field's logical volume, and return the transform of that placement. +G4ThreeVector GateField::findContainingPlacement( + const G4ThreeVector &worldPoint, + const G4AffineTransform *&outTransform +) const { + + // Loop over all placements once: + // - return immediately if the point is inside any placement; + // - otherwise accumulate the closest surface distance for the fallback. + // This avoids re-transforming the point in a second pass. + std::size_t closestIdx = 0; + double minDistToSurface = std::numeric_limits::infinity(); + G4ThreeVector closestLocal{}; + + for (std::size_t i = 0; i < m_transforms.size(); ++i) { + const auto& tr = m_transforms[i]; + const G4ThreeVector localPoint = tr.InverseTransformPoint(worldPoint); + + if (m_solid->Inside(localPoint) != kOutside) { + outTransform = &tr; + return localPoint; + } + + // Fallback: track the placement whose surface is nearest. + // DistanceToIn is valid here because Inside() just returned kOutside. + const double d = m_solid->DistanceToIn(localPoint); + if (d < minDistToSurface) { + minDistToSurface = d; + closestIdx = i; + closestLocal = localPoint; + } + } + + // sanity check: if the closest surface is still too far, this is likely a real bug + if (minDistToSurface > m_fallbackFatalDistanceMM) { + std::ostringstream msg; + msg << "GateField::findContainingPlacement: world point (" << worldPoint.x() + << ", " << worldPoint.y() << ", " << worldPoint.z() << ") mm is " + << minDistToSurface + << " mm outside every cached placement of the field's solid — " + << "well beyond any chord-finder overshoot.\n" + << " Closest placement: index " << closestIdx << " (local point " + << closestLocal.x() << ", " << closestLocal.y() << ", " + << closestLocal.z() << ").\n" + << " Maximum allowed distance before fatal: " << m_fallbackFatalDistanceMM << " mm.\n" + << " This likely indicates a real bug in the geometry or field setup\n"; + + G4Exception( + "GateField::findContainingPlacement", + "GateField0001", + FatalException, + msg.str().c_str() + ); + } + + outTransform = &m_transforms[closestIdx]; + + return closestLocal; + +} + +// rotate a field vector from local to world coordinates using the given +// transform +G4ThreeVector GateField::rotateToWorld(const G4ThreeVector &localField, + const G4AffineTransform &transform) { + return transform.TransformAxis(localField); +} diff --git a/core/opengate_core/opengate_lib/GateField.h b/core/opengate_core/opengate_lib/GateField.h new file mode 100644 index 0000000000..cee3bd7471 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateField.h @@ -0,0 +1,69 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateField_h +#define GateField_h + +#include "G4AffineTransform.hh" +#include "G4RotationMatrix.hh" +#include "G4ThreeVector.hh" +#include + +// Forward declaration to avoid including the full header. +class G4VSolid; + +// Helper base class for fields in GATE: +// - stores the world-to-local transforms of every physical +// placement of the logical volume the field is attached to. +// - provides the coordinate conversions needed to evaluate fields in each +// placement's local frame and rotate the result back to world. +class GateField { + +public: + // constructor — deltaChordMM is the chord-finder tolerance (delta_chord), + // used to compute the fallback-fatal distance internally. + GateField( + const G4VSolid *solid, + std::vector translations, + std::vector rotations, + double deltaChordMM + ); + + // update the transforms (e.g. after a geometry change between runs) + void SetTransforms( + std::vector translations, + std::vector rotations + ); + +protected: + // find the local coordinates of worldPoint in the containing placement, and + // return the world-to-local transform of that placement in outTransform. + G4ThreeVector findContainingPlacement( + const G4ThreeVector &worldPoint, + const G4AffineTransform *&outTransform + ) const; + + // rotate a field vector from local to world coordinates using the given + // transform. + static G4ThreeVector rotateToWorld( + const G4ThreeVector &localField, + const G4AffineTransform &transform + ); + + const G4VSolid *m_solid; // solid of the logical volume the field is attached to + + double m_fallbackFatalDistanceMM; // computed from deltaChordMM at construction + + std::vector m_transforms; // one per physical placement + + // Upper bound on the radius of curvature [mm] used to derive + // the fallback-fatal distance. + static constexpr double kMaxCurvatureRadiusMM = 100.0e3; + +}; + +#endif // GateField_h diff --git a/core/opengate_core/opengate_lib/GateMagneticField.cpp b/core/opengate_core/opengate_lib/GateMagneticField.cpp new file mode 100644 index 0000000000..947a2e2a6d --- /dev/null +++ b/core/opengate_core/opengate_lib/GateMagneticField.cpp @@ -0,0 +1,50 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include "GateMagneticField.h" + +// constructor +GateMagneticField::GateMagneticField( + G4MagneticField *inner, + const G4VSolid *solid, + std::vector translations, + std::vector rotations, + double deltaChordMM +) : G4MagneticField(), GateField(solid, std::move(translations), std::move(rotations), deltaChordMM), m_inner(inner) +{ + +} + +// override of G4MagneticField interface +void GateMagneticField::GetFieldValue( + const G4double Point[4], + G4double *Bfield +) const +{ + + // get the global coordinates of the point and the transform of the containing placement + const G4ThreeVector worldPoint(Point[0], Point[1], Point[2]); + + // get the local coordinates of the point in the containing placement, and the transform of that placement + const G4AffineTransform *transform = nullptr; + const G4ThreeVector localPoint = + findContainingPlacement(worldPoint, transform); + const G4double localPos[4] = {localPoint.x(), localPoint.y(), localPoint.z(), Point[3]}; + + // get the local field value from the inner field object + G4double localB[3] = {0.0, 0.0, 0.0}; + m_inner->GetFieldValue(localPos, localB); + + // rotate the field back to world coordinates + const G4ThreeVector worldB = rotateToWorld({localB[0], localB[1], localB[2]}, *transform); + + // copy the result to the output array + Bfield[0] = worldB.x(); + Bfield[1] = worldB.y(); + Bfield[2] = worldB.z(); + +} diff --git a/core/opengate_core/opengate_lib/GateMagneticField.h b/core/opengate_core/opengate_lib/GateMagneticField.h new file mode 100644 index 0000000000..d1c3a230d6 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateMagneticField.h @@ -0,0 +1,58 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateMagneticField_h +#define GateMagneticField_h + +#include "G4MagneticField.hh" +#include "G4RotationMatrix.hh" +#include "G4ThreeVector.hh" +#include "GateField.h" +#include + +class G4VSolid; + +// GATE wrapper for G4MagneticField that transforms the query point to the local +// coordinates of the physical volume(s) the field is attached to, then +// delegates to the inner field class to get the field value, and finally +// transforms the field vector back to world coordinates. + +// wrapper for an internal G4MagneticField converted to local coordinates +class GateMagneticField : public G4MagneticField, protected GateField +{ +public: + // constructor — deltaChordMM is the chord-finder tolerance (delta_chord), + // forwarded to GateField to compute the fallback-fatal distance internally. + GateMagneticField( + G4MagneticField *inner, // wrapped inner field + const G4VSolid *solid, // solid of the logical volume the field is attached to + std::vector translations, // translations of the physical placements + std::vector rotations, // rotations of the physical placements of the logical volume + double deltaChordMM + ); + + // override of G4MagneticField interface + void GetFieldValue( + const G4double Point[4], + G4double *Bfield + ) const override; + + // forward override of GateField::SetTransforms + inline void SetTransforms( + std::vector translations, + std::vector rotations + ) + { + GateField::SetTransforms(std::move(translations), std::move(rotations)); + } + +private: + G4MagneticField *m_inner; // wrapped inner field + +}; + +#endif // GateMagneticField_h diff --git a/core/opengate_core/opengate_lib/pyGateMagneticField.cpp b/core/opengate_core/opengate_lib/pyGateMagneticField.cpp new file mode 100644 index 0000000000..da23ddd6a4 --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateMagneticField.cpp @@ -0,0 +1,39 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include +#include + +namespace py = pybind11; + +#include "G4MagneticField.hh" +#include "G4VSolid.hh" +#include "GateMagneticField.h" + +// python bindings for GateMagneticField +void init_GateMagneticField(py::module &m) +{ + + py::class_>( + m, "GateMagneticField") + + .def(py::init([](G4MagneticField *inner, const G4VSolid *solid, + std::vector translations, + std::vector rotations, + double delta_chord_mm) { + return new GateMagneticField(inner, solid, translations, + rotations, delta_chord_mm); + }), + py::arg("inner_field"), py::arg("solid"), py::arg("translations"), + py::arg("rotations"), py::arg("delta_chord_mm")) + + .def("SetTransforms", &GateMagneticField::SetTransforms, + py::arg("translations"), py::arg("rotations"), + "Replace the cached world-to-local transforms. Call between runs " + "only after dynamic geometry changes have been applied."); +} From 64fe1ba8d171b8e48e936cc17060f6187a8504fa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Tue, 5 May 2026 10:46:22 +0200 Subject: [PATCH 02/28] Enhance field management on the Python side by adding support for field definition in local coordinates and updates for dynamic geometry changes --- opengate/actors/dynamicactors.py | 17 +++ opengate/engines.py | 1 + opengate/geometry/fields.py | 245 +++++++++++++++++++++++++++---- 3 files changed, 235 insertions(+), 28 deletions(-) diff --git a/opengate/actors/dynamicactors.py b/opengate/actors/dynamicactors.py index cf89573c6c..07112fcdad 100644 --- a/opengate/actors/dynamicactors.py +++ b/opengate/actors/dynamicactors.py @@ -60,6 +60,23 @@ def BeginOfRunActionMasterThread(self, run_id): # CloseGeometry: pOptimise=true, verbose=false, G4VPhysicalVolume *vol=0 gm.CloseGeometry(self.simulation.dyn_geom_optimise, False, None) + # Refresh field transforms after geometry changes + # This ensures fields attached to moving volumes use the updated transforms + self._refresh_field_transforms(run_id) + + def _refresh_field_transforms(self, run_id): + """Recompute world-to-local transforms for every registered field. + + Called after the geometry-changers for `run_id` have been applied, + so that fields attached to moving volumes see the updated placement + transforms in the upcoming run. + """ + fields = getattr(self.simulation.volume_manager, "fields", None) + if not fields: + return + for field in fields.values(): + field.refresh_transforms() + class DynamicSourceActor(DynamicActorBase, g4.GateVActor): diff --git a/opengate/engines.py b/opengate/engines.py index 8ff27c26d2..a5f8b7b0fd 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -953,6 +953,7 @@ def ConstructSDandField(self): for field in self.simulation_engine.simulation.volume_manager.fields.values(): for volume_name in field.attached_to: volume_obj = self.volume_manager.get_volume(volume_name) + field._field_volume_obj = volume_obj volume_obj.g4_field_manager = field.create_field_manager() volume_obj.g4_logical_volume.SetFieldManager( volume_obj.g4_field_manager, True diff --git a/opengate/geometry/fields.py b/opengate/geometry/fields.py index a487bcb91f..506de526d9 100644 --- a/opengate/geometry/fields.py +++ b/opengate/geometry/fields.py @@ -1,12 +1,21 @@ +from typing import Any + import opengate_core as g4 from ..base import GateObject, process_cls +from ..geometry.utility import ( + get_transform_world_to_local, + vec_np_as_g4, + rot_np_as_g4, + vec_g4_as_np, + rot_g4_as_np, +) from ..utility import g4_units # ! ======= KNOWN TODO'S ======== # ! - implement the possibility of choosing the stepper type and equation type # ! - bind the sextupole magnetic field geant4 implementation -# ! - implement mapped fields (e.g., from a CSV file) +# ! WIP - implement mapped fields (e.g., from a CSV file) # ! - Overhead for custom fields implementation: every GetFieldValue call crosses c++ -> Python -> c++, acquiring the GIL each time, # ! which is very inefficient. Need to implement a more efficient way on the C++ side. # ! - in MT mode, create_field_manager() is called per thread and overwrites shared instance attrs (g4_field, etc.), @@ -70,6 +79,7 @@ def __init__(self, *args, **kwargs) -> None: self.g4_field = None self._g4_runtime_objects = [] + self._field_volume_obj: Any = None # set by engines.py before create_field_manager() is called self.attached_to = [] self._field_changes_energy = False @@ -79,6 +89,33 @@ def field_changes_energy(self) -> bool: """Whether the field changes particle energy (False for magnetic, True for others).""" return self._field_changes_energy + def refresh_transforms(self) -> None: + """Recompute and push cached world-to-local transforms after dynamic + geometry changes. + """ + if self.g4_field is None or self._field_volume_obj is None: + return + + volume = self._field_volume_obj + g4_translations = [] + g4_rotations = [] + for i in range(volume.number_of_repetitions): + pv = volume.get_g4_physical_volume(i) + T = vec_g4_as_np(pv.GetObjectTranslation()) + R = rot_g4_as_np(pv.GetObjectRotation()) + for anc in volume.ancestor_volumes[::-1]: + anc_pv = getattr(anc, "g4_physical_volume", None) + if anc_pv is None: + continue + anc_T = vec_g4_as_np(anc_pv.GetObjectTranslation()) + anc_R = rot_g4_as_np(anc_pv.GetObjectRotation()) + T = anc_R @ T + anc_T + R = anc_R @ R + g4_translations.append(vec_np_as_g4(T)) + g4_rotations.append(rot_np_as_g4(R)) + + self.g4_field.SetTransforms(g4_translations, g4_rotations) + def create_field_manager(self) -> g4.G4FieldManager: """Construct the field and return a configured G4FieldManager.""" msg = "create_field_manager() must be implemented in subclasses." @@ -104,40 +141,50 @@ class MagneticField(FieldBase): def __init__(self, *args, **kwargs) -> None: super().__init__(*args, **kwargs) - # Magnetic fields don't change particle energy self._field_changes_energy = False - self.g4_equation_of_motion = None self.g4_integrator_stepper = None self.g4_integration_driver = None self.g4_chord_finder = None - def _create_field(self) -> None: - """Create the G4 field object. Override.""" - msg = "_create_field() must be implemented in subclasses." - raise NotImplementedError(msg) + def _create_inner_field(self): + """Create and return the inner G4MagneticField in local volume coordinates. + Override in subclasses.""" + raise NotImplementedError("_create_inner_field() must be implemented in subclasses.") def create_field_manager(self) -> g4.G4FieldManager: """Construct the field and return a configured G4FieldManager.""" - # Create the G4 field object (subclass responsibility) - self._create_field() + # Build the inner field (subclass responsibility, in local volume coords). + inner = self._create_inner_field() - # Create equation of motion, stepper, driver, chord finder TODO: allow user choice + # Wrap in GateMagneticField: handles multi-placement coordinate transforms + # and the integration-overshoot fallback uniformly for all field types. + translations_np, rotations_np = get_transform_world_to_local( + self._field_volume_obj + ) + self.g4_field = g4.GateMagneticField( + inner, + self._field_volume_obj.g4_solid, + [vec_np_as_g4(t) for t in translations_np], + [rot_np_as_g4(r) for r in rotations_np], + self.delta_chord, + ) + + # TODO: allow user choice of stepper and equation type self.g4_equation_of_motion = g4.G4Mag_UsualEqRhs(self.g4_field) self.g4_integrator_stepper = g4.G4ClassicalRK4( self.g4_equation_of_motion, - 6, # number of variables for magnetic field = 6 (x,y,z + px,py,pz) + 6, # x,y,z + px,py,pz ) self.g4_integration_driver = g4.G4MagInt_Driver( self.step_minimum, self.g4_integrator_stepper, - 6, # number of variables for magnetic field = 6 (x,y,z + px,py,pz) + 6, 0, # no verbosity ) self.g4_chord_finder = g4.G4ChordFinder(self.g4_integration_driver) self.g4_chord_finder.SetDeltaChord(self.delta_chord) - # Create and configure field manager fm = g4.G4FieldManager( self.g4_field, self.g4_chord_finder, self.field_changes_energy ) @@ -146,9 +193,12 @@ def create_field_manager(self) -> g4.G4FieldManager: fm.SetMinimumEpsilonStep(self.min_epsilon_step) fm.SetMaximumEpsilonStep(self.max_epsilon_step) - # Keep runtime objects alive for the full simulation lifetime. + # Keep all runtime objects alive for the full simulation lifetime. + # 'inner_field' is held as a raw pointer by GateMagneticField in C++, + # so Python must retain a reference here. self._g4_runtime_objects.append( { + "inner_field": inner, "field": self.g4_field, "equation": self.g4_equation_of_motion, "stepper": self.g4_integrator_stepper, @@ -169,7 +219,12 @@ def close(self) -> None: class UniformMagneticField(MagneticField): - """Uniform magnetic field with constant field vector.""" + """Uniform magnetic field with constant field vector. + + field_vector is specified in the local coordinate frame of the attached + volume. For a non-rotated volume this is identical to the world frame. + For a rotated volume the field direction rotates with the volume. + """ # hints for IDE field_vector: list @@ -178,7 +233,8 @@ class UniformMagneticField(MagneticField): "field_vector": ( [0, 0, 0], { - "doc": "Field vector [Bx, By, Bz]. Each component in magnetic field strength units (e.g., Tesla).", + "doc": "Field vector [Bx, By, Bz] in local volume coordinates. " + "Each component in magnetic field strength units (e.g., Tesla).", }, ), } @@ -186,9 +242,8 @@ class UniformMagneticField(MagneticField): def __init__(self, *args, **kwargs) -> None: super().__init__(*args, **kwargs) - def _create_field(self) -> None: - """Create the uniform magnetic field.""" - self.g4_field = g4.G4UniformMagField(g4.G4ThreeVector(*self.field_vector)) + def _create_inner_field(self): + return g4.G4UniformMagField(g4.G4ThreeVector(*self.field_vector)) class QuadrupoleMagneticField(MagneticField): @@ -209,9 +264,8 @@ class QuadrupoleMagneticField(MagneticField): def __init__(self, *args, **kwargs) -> None: super().__init__(*args, **kwargs) - def _create_field(self) -> None: - """Create the quadrupole magnetic field.""" - self.g4_field = g4.G4QuadrupoleMagField(self.gradient) + def _create_inner_field(self): + return g4.G4QuadrupoleMagField(self.gradient) class CustomMagneticField(MagneticField): @@ -232,22 +286,25 @@ class CustomMagneticField(MagneticField): def __init__(self, *args, **kwargs) -> None: super().__init__(*args, **kwargs) - def _create_field(self) -> None: - """Create the custom magnetic field using the Python trampoline.""" + def _create_inner_field(self): + """Create the custom magnetic field using the Python trampoline. + + field_function receives (x, y, z, t) in the local coordinate frame of + the attached volume and must return [Bx, By, Bz] in that same local + frame. The base class rotates the result to world coordinates. + """ if self.field_function is None: raise ValueError("field_function must be provided for CustomMagneticField") if not callable(self.field_function): raise TypeError("field_function must be a callable function") - # Check if the function returns 3 components test_point = [0, 0, 0, 0] - result = self.field_function(*test_point) + result = list(self.field_function(*test_point)) if len(result) != 3: raise ValueError( "field_function must return a list of 3 components: [Bx, By, Bz]" ) - # Create a custom G4MagneticField subclass that calls our Python function class _PyMagneticField(g4.G4MagneticField): def __init__(inner_self, callback): super().__init__() @@ -256,7 +313,7 @@ def __init__(inner_self, callback): def GetFieldValue(inner_self, point): return inner_self._callback(*point) - self.g4_field = _PyMagneticField(self.field_function) + return _PyMagneticField(self.field_function) def to_dictionary(self): raise NotImplementedError( @@ -530,10 +587,141 @@ def to_dictionary(self): ) +# class MappedMagneticField(MagneticField): +# """Magnetic field defined by values on a regular 3D Cartesian grid. + +# Field lookup is performed entirely in C++ (no Python callbacks at tracking +# time) using either trilinear or nearest-neighbour interpolation. Points +# outside the grid extent return zero field. +# """ + +# # ! TODO's (@srmarcballestero): +# # ! - refactor the _create_field() function to make it simpler. +# # ! - warning when using nearest interpolation with coarse grids +# # ! - implement a MappedElectroMagneticField, and MappedElectricField. +# # ! + for the electromagnetic case, separate grids should be used. +# # ! + a mother class should abstract the core functionality. + +# # hints for IDE +# field_matrix: np.ndarray +# interpolation: str + +# user_info_defaults = { +# "field_matrix": ( +# None, +# { +# "doc": ( +# "2D numpy array of shape (N, 6) with columns [x, y, z, Bx, By, Bz] " +# "on a regular Cartesian grid. Geant4 units should be used." +# ), +# }, +# ), +# "interpolation": ( +# "trilinear", +# { +# "doc": "Interpolation method: 'trilinear' (default) or 'nearest'.", +# }, +# ), +# } + +# def __init__(self, *args, **kwargs) -> None: +# super().__init__(*args, **kwargs) + +# def _create_field(self) -> None: +# if self.field_matrix is None: +# raise ValueError("field_matrix must be provided for MappedMagneticField") + +# mat = np.asarray(self.field_matrix, dtype=np.float64) +# if mat.ndim != 2 or mat.shape[1] != 6: +# raise ValueError( +# "field_matrix must be a 2D array with shape (N, 6) " +# "containing columns [x, y, z, Bx, By, Bz]" +# ) + +# # Separate coordinates and field values +# positions = mat[:, :3] +# B_values = mat[:, 3:] + +# # Sort in lexicographical order: x slowest, z fastest — matches C++ flat index. +# sort_idx = np.lexsort((positions[:, 2], positions[:, 1], positions[:, 0])) +# positions = positions[sort_idx] +# B_values = B_values[sort_idx] + +# # Round to suppress floating-point noise before uniqueness check. +# x_vals = np.unique(np.round(positions[:, 0], 10)) +# y_vals = np.unique(np.round(positions[:, 1], 10)) +# z_vals = np.unique(np.round(positions[:, 2], 10)) +# nx, ny, nz = len(x_vals), len(y_vals), len(z_vals) + +# # --- Grid validation --- +# if len(positions) != nx * ny * nz: +# raise ValueError( +# f"field_matrix does not define a complete regular 3D grid: " +# f"expected {nx}*{ny}*{nz}={nx * ny * nz} points, got {len(positions)}" +# ) +# if nx > 2 and not np.allclose(np.diff(x_vals), x_vals[1] - x_vals[0]): +# raise ValueError("field_matrix x-axis does not have uniform spacing") +# if ny > 2 and not np.allclose(np.diff(y_vals), y_vals[1] - y_vals[0]): +# raise ValueError("field_matrix y-axis does not have uniform spacing") +# if nz > 2 and not np.allclose(np.diff(z_vals), z_vals[1] - z_vals[0]): +# raise ValueError("field_matrix z-axis does not have uniform spacing") +# if nx < 2 or ny < 2 or nz < 2: +# raise ValueError( +# "field_matrix must have at least 2 unique points along each axis " +# "(no degenerate axes allowed)" +# ) +# # ---------------------- + +# x0, y0, z0 = float(x_vals[0]), float(y_vals[0]), float(z_vals[0]) +# dx = float(x_vals[1] - x_vals[0]) +# dy = float(y_vals[1] - y_vals[0]) +# dz = float(z_vals[1] - z_vals[0]) + +# interp_map = { +# "trilinear": g4.GateMappedMagneticFieldInterpolation.Trilinear, +# "nearest": g4.GateMappedMagneticFieldInterpolation.Nearest, +# } +# if self.interpolation not in interp_map: +# raise ValueError( +# f"Unknown interpolation '{self.interpolation}'. " +# f"Available options are 'trilinear' or 'nearest'." +# ) + +# # Collect one local-to-world transform per physical placement. +# translations_np, rotations_np = get_transform_world_to_local( +# self._field_volume_obj +# ) + +# g4_translations = [vec_np_as_g4(t) for t in translations_np] +# g4_rotations = [rot_np_as_g4(r) for r in rotations_np] + +# self.g4_field = g4.GateMappedMagneticField( +# B_values[:, 0], B_values[:, 1], B_values[:, 2], +# nx, ny, nz, +# x0, y0, z0, +# dx, dy, dz, +# interp_map[self.interpolation], +# g4_translations, g4_rotations, +# ) + +# # Serialization +# # def to_dictionary(self): +# # d = super().to_dictionary() +# # if self.field_matrix is not None: +# # d["field_matrix"] = np.asarray(self.field_matrix).tolist() +# # return d + +# # def from_dictionary(self, d): +# # super().from_dictionary(d) +# # if "field_matrix" in d and d["field_matrix"] is not None: +# # self.field_matrix = np.asarray(d["field_matrix"]) + + field_types = { "UniformMagneticField": UniformMagneticField, "QuadrupoleMagneticField": QuadrupoleMagneticField, "CustomMagneticField": CustomMagneticField, + # "MappedMagneticField": MappedMagneticField, "UniformElectricField": UniformElectricField, "CustomElectricField": CustomElectricField, "UniformElectroMagneticField": UniformElectroMagneticField, @@ -545,6 +733,7 @@ def to_dictionary(self): process_cls(UniformMagneticField) process_cls(QuadrupoleMagneticField) process_cls(CustomMagneticField) +# process_cls(MappedMagneticField) process_cls(ElectroMagneticField) process_cls(ElectricField) process_cls(UniformElectricField) From 7e3a60191dd6ad9d035dd443941172b34b514877 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 5 May 2026 08:55:12 +0000 Subject: [PATCH 03/28] [pre-commit.ci] Automatic python and c++ formatting --- core/opengate_core/opengate_lib/GateField.cpp | 55 ++++++++----------- core/opengate_core/opengate_lib/GateField.h | 34 +++++------- .../opengate_lib/GateMagneticField.cpp | 38 ++++++------- .../opengate_lib/GateMagneticField.h | 32 +++++------ .../opengate_lib/pyGateMagneticField.cpp | 7 +-- opengate/geometry/fields.py | 8 ++- 6 files changed, 77 insertions(+), 97 deletions(-) diff --git a/core/opengate_core/opengate_lib/GateField.cpp b/core/opengate_core/opengate_lib/GateField.cpp index eb8367b8d9..7a1c4ec572 100644 --- a/core/opengate_core/opengate_lib/GateField.cpp +++ b/core/opengate_core/opengate_lib/GateField.cpp @@ -11,21 +11,20 @@ #include "globals.hh" // G4Exception #include +#include #include #include #include -#include #include // constructor -GateField::GateField( - const G4VSolid *solid, - std::vector translations, - std::vector rotations, - double deltaChordMM -) : m_solid(solid), - m_fallbackFatalDistanceMM(5.0 * std::sqrt(8.0 * kMaxCurvatureRadiusMM * deltaChordMM)) -{ +GateField::GateField(const G4VSolid *solid, + std::vector translations, + std::vector rotations, + double deltaChordMM) + : m_solid(solid), + m_fallbackFatalDistanceMM( + 5.0 * std::sqrt(8.0 * kMaxCurvatureRadiusMM * deltaChordMM)) { // sanity-check the inputs before caching if (solid == nullptr) throw std::invalid_argument("GateField: solid must not be null"); @@ -34,19 +33,17 @@ GateField::GateField( throw std::invalid_argument("GateField: translations and rotations must be " "non-empty and have the same size"); - - // initial cache the world-to-local transforms for every physical placement of the physical volume + // initial cache the world-to-local transforms for every physical placement of + // the physical volume m_transforms.reserve(translations.size()); for (std::size_t i = 0; i < translations.size(); ++i) m_transforms.emplace_back(rotations[i].inverse(), translations[i]); - } -// recache the world-to-local transforms (e.g. after a geometry change between runs) -void GateField::SetTransforms( - std::vector translations, - std::vector rotations -) { +// recache the world-to-local transforms (e.g. after a geometry change between +// runs) +void GateField::SetTransforms(std::vector translations, + std::vector rotations) { // sanity-check the inputs before caching if (translations.size() != rotations.size() || translations.empty()) @@ -54,20 +51,19 @@ void GateField::SetTransforms( "GateField::SetTransforms: translations and rotations must be " "non-empty and have the same size"); - // recache the world-to-local transforms for every physical placement of the logical volume + // recache the world-to-local transforms for every physical placement of the + // logical volume m_transforms.clear(); m_transforms.reserve(translations.size()); for (std::size_t i = 0; i < translations.size(); ++i) m_transforms.emplace_back(rotations[i].inverse(), translations[i]); - } // find the local coordinates of worldPoint in the containing placement of the // field's logical volume, and return the transform of that placement. G4ThreeVector GateField::findContainingPlacement( const G4ThreeVector &worldPoint, - const G4AffineTransform *&outTransform -) const { + const G4AffineTransform *&outTransform) const { // Loop over all placements once: // - return immediately if the point is inside any placement; @@ -78,7 +74,7 @@ G4ThreeVector GateField::findContainingPlacement( G4ThreeVector closestLocal{}; for (std::size_t i = 0; i < m_transforms.size(); ++i) { - const auto& tr = m_transforms[i]; + const auto &tr = m_transforms[i]; const G4ThreeVector localPoint = tr.InverseTransformPoint(worldPoint); if (m_solid->Inside(localPoint) != kOutside) { @@ -96,7 +92,8 @@ G4ThreeVector GateField::findContainingPlacement( } } - // sanity check: if the closest surface is still too far, this is likely a real bug + // sanity check: if the closest surface is still too far, this is likely a + // real bug if (minDistToSurface > m_fallbackFatalDistanceMM) { std::ostringstream msg; msg << "GateField::findContainingPlacement: world point (" << worldPoint.x() @@ -107,21 +104,17 @@ G4ThreeVector GateField::findContainingPlacement( << " Closest placement: index " << closestIdx << " (local point " << closestLocal.x() << ", " << closestLocal.y() << ", " << closestLocal.z() << ").\n" - << " Maximum allowed distance before fatal: " << m_fallbackFatalDistanceMM << " mm.\n" + << " Maximum allowed distance before fatal: " + << m_fallbackFatalDistanceMM << " mm.\n" << " This likely indicates a real bug in the geometry or field setup\n"; - G4Exception( - "GateField::findContainingPlacement", - "GateField0001", - FatalException, - msg.str().c_str() - ); + G4Exception("GateField::findContainingPlacement", "GateField0001", + FatalException, msg.str().c_str()); } outTransform = &m_transforms[closestIdx]; return closestLocal; - } // rotate a field vector from local to world coordinates using the given diff --git a/core/opengate_core/opengate_lib/GateField.h b/core/opengate_core/opengate_lib/GateField.h index cee3bd7471..66b002c0e1 100644 --- a/core/opengate_core/opengate_lib/GateField.h +++ b/core/opengate_core/opengate_lib/GateField.h @@ -26,44 +26,36 @@ class GateField { public: // constructor — deltaChordMM is the chord-finder tolerance (delta_chord), // used to compute the fallback-fatal distance internally. - GateField( - const G4VSolid *solid, - std::vector translations, - std::vector rotations, - double deltaChordMM - ); + GateField(const G4VSolid *solid, std::vector translations, + std::vector rotations, double deltaChordMM); // update the transforms (e.g. after a geometry change between runs) - void SetTransforms( - std::vector translations, - std::vector rotations - ); + void SetTransforms(std::vector translations, + std::vector rotations); protected: // find the local coordinates of worldPoint in the containing placement, and // return the world-to-local transform of that placement in outTransform. - G4ThreeVector findContainingPlacement( - const G4ThreeVector &worldPoint, - const G4AffineTransform *&outTransform - ) const; + G4ThreeVector + findContainingPlacement(const G4ThreeVector &worldPoint, + const G4AffineTransform *&outTransform) const; // rotate a field vector from local to world coordinates using the given // transform. - static G4ThreeVector rotateToWorld( - const G4ThreeVector &localField, - const G4AffineTransform &transform - ); + static G4ThreeVector rotateToWorld(const G4ThreeVector &localField, + const G4AffineTransform &transform); - const G4VSolid *m_solid; // solid of the logical volume the field is attached to + const G4VSolid + *m_solid; // solid of the logical volume the field is attached to - double m_fallbackFatalDistanceMM; // computed from deltaChordMM at construction + double + m_fallbackFatalDistanceMM; // computed from deltaChordMM at construction std::vector m_transforms; // one per physical placement // Upper bound on the radius of curvature [mm] used to derive // the fallback-fatal distance. static constexpr double kMaxCurvatureRadiusMM = 100.0e3; - }; #endif // GateField_h diff --git a/core/opengate_core/opengate_lib/GateMagneticField.cpp b/core/opengate_core/opengate_lib/GateMagneticField.cpp index 947a2e2a6d..d35f0035a0 100644 --- a/core/opengate_core/opengate_lib/GateMagneticField.cpp +++ b/core/opengate_core/opengate_lib/GateMagneticField.cpp @@ -8,43 +8,41 @@ #include "GateMagneticField.h" // constructor -GateMagneticField::GateMagneticField( - G4MagneticField *inner, - const G4VSolid *solid, - std::vector translations, - std::vector rotations, - double deltaChordMM -) : G4MagneticField(), GateField(solid, std::move(translations), std::move(rotations), deltaChordMM), m_inner(inner) -{ - -} +GateMagneticField::GateMagneticField(G4MagneticField *inner, + const G4VSolid *solid, + std::vector translations, + std::vector rotations, + double deltaChordMM) + : G4MagneticField(), GateField(solid, std::move(translations), + std::move(rotations), deltaChordMM), + m_inner(inner) {} // override of G4MagneticField interface -void GateMagneticField::GetFieldValue( - const G4double Point[4], - G4double *Bfield -) const -{ +void GateMagneticField::GetFieldValue(const G4double Point[4], + G4double *Bfield) const { - // get the global coordinates of the point and the transform of the containing placement + // get the global coordinates of the point and the transform of the containing + // placement const G4ThreeVector worldPoint(Point[0], Point[1], Point[2]); - // get the local coordinates of the point in the containing placement, and the transform of that placement + // get the local coordinates of the point in the containing placement, and the + // transform of that placement const G4AffineTransform *transform = nullptr; const G4ThreeVector localPoint = findContainingPlacement(worldPoint, transform); - const G4double localPos[4] = {localPoint.x(), localPoint.y(), localPoint.z(), Point[3]}; + const G4double localPos[4] = {localPoint.x(), localPoint.y(), localPoint.z(), + Point[3]}; // get the local field value from the inner field object G4double localB[3] = {0.0, 0.0, 0.0}; m_inner->GetFieldValue(localPos, localB); // rotate the field back to world coordinates - const G4ThreeVector worldB = rotateToWorld({localB[0], localB[1], localB[2]}, *transform); + const G4ThreeVector worldB = + rotateToWorld({localB[0], localB[1], localB[2]}, *transform); // copy the result to the output array Bfield[0] = worldB.x(); Bfield[1] = worldB.y(); Bfield[2] = worldB.z(); - } diff --git a/core/opengate_core/opengate_lib/GateMagneticField.h b/core/opengate_core/opengate_lib/GateMagneticField.h index d1c3a230d6..53bff33896 100644 --- a/core/opengate_core/opengate_lib/GateMagneticField.h +++ b/core/opengate_core/opengate_lib/GateMagneticField.h @@ -22,37 +22,31 @@ class G4VSolid; // transforms the field vector back to world coordinates. // wrapper for an internal G4MagneticField converted to local coordinates -class GateMagneticField : public G4MagneticField, protected GateField -{ +class GateMagneticField : public G4MagneticField, protected GateField { public: // constructor — deltaChordMM is the chord-finder tolerance (delta_chord), // forwarded to GateField to compute the fallback-fatal distance internally. - GateMagneticField( - G4MagneticField *inner, // wrapped inner field - const G4VSolid *solid, // solid of the logical volume the field is attached to - std::vector translations, // translations of the physical placements - std::vector rotations, // rotations of the physical placements of the logical volume - double deltaChordMM - ); + GateMagneticField(G4MagneticField *inner, // wrapped inner field + const G4VSolid *solid, // solid of the logical volume the + // field is attached to + std::vector + translations, // translations of the physical placements + std::vector + rotations, // rotations of the physical placements of + // the logical volume + double deltaChordMM); // override of G4MagneticField interface - void GetFieldValue( - const G4double Point[4], - G4double *Bfield - ) const override; + void GetFieldValue(const G4double Point[4], G4double *Bfield) const override; // forward override of GateField::SetTransforms - inline void SetTransforms( - std::vector translations, - std::vector rotations - ) - { + inline void SetTransforms(std::vector translations, + std::vector rotations) { GateField::SetTransforms(std::move(translations), std::move(rotations)); } private: G4MagneticField *m_inner; // wrapped inner field - }; #endif // GateMagneticField_h diff --git a/core/opengate_core/opengate_lib/pyGateMagneticField.cpp b/core/opengate_core/opengate_lib/pyGateMagneticField.cpp index da23ddd6a4..c1a26b588c 100644 --- a/core/opengate_core/opengate_lib/pyGateMagneticField.cpp +++ b/core/opengate_core/opengate_lib/pyGateMagneticField.cpp @@ -15,8 +15,7 @@ namespace py = pybind11; #include "GateMagneticField.h" // python bindings for GateMagneticField -void init_GateMagneticField(py::module &m) -{ +void init_GateMagneticField(py::module &m) { py::class_>( @@ -26,8 +25,8 @@ void init_GateMagneticField(py::module &m) std::vector translations, std::vector rotations, double delta_chord_mm) { - return new GateMagneticField(inner, solid, translations, - rotations, delta_chord_mm); + return new GateMagneticField(inner, solid, translations, rotations, + delta_chord_mm); }), py::arg("inner_field"), py::arg("solid"), py::arg("translations"), py::arg("rotations"), py::arg("delta_chord_mm")) diff --git a/opengate/geometry/fields.py b/opengate/geometry/fields.py index 506de526d9..0b7a6b67a4 100644 --- a/opengate/geometry/fields.py +++ b/opengate/geometry/fields.py @@ -79,7 +79,9 @@ def __init__(self, *args, **kwargs) -> None: self.g4_field = None self._g4_runtime_objects = [] - self._field_volume_obj: Any = None # set by engines.py before create_field_manager() is called + self._field_volume_obj: Any = ( + None # set by engines.py before create_field_manager() is called + ) self.attached_to = [] self._field_changes_energy = False @@ -150,7 +152,9 @@ def __init__(self, *args, **kwargs) -> None: def _create_inner_field(self): """Create and return the inner G4MagneticField in local volume coordinates. Override in subclasses.""" - raise NotImplementedError("_create_inner_field() must be implemented in subclasses.") + raise NotImplementedError( + "_create_inner_field() must be implemented in subclasses." + ) def create_field_manager(self) -> g4.G4FieldManager: """Construct the field and return a configured G4FieldManager.""" From 2a552eaad6782d0427621d980140df2f8347777d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Tue, 5 May 2026 17:37:43 +0200 Subject: [PATCH 04/28] Enhance support for electromagnetic fields (E or E+B) , and implement a faster UniformElectroMagneticField --- core/opengate_core/opengate_core.cpp | 4 + .../opengate_lib/GateElectroMagneticField.cpp | 42 ++ .../opengate_lib/GateElectroMagneticField.h | 41 ++ .../GateUniformElectroMagneticField.cpp | 24 ++ .../GateUniformElectroMagneticField.h | 32 ++ .../pyGateElectroMagneticField.cpp | 36 ++ .../pyGateUniformElectroMagneticField.cpp | 27 ++ opengate/engines.py | 1 + opengate/geometry/fields.py | 372 ++++++------------ 9 files changed, 338 insertions(+), 241 deletions(-) create mode 100644 core/opengate_core/opengate_lib/GateElectroMagneticField.cpp create mode 100644 core/opengate_core/opengate_lib/GateElectroMagneticField.h create mode 100644 core/opengate_core/opengate_lib/GateUniformElectroMagneticField.cpp create mode 100644 core/opengate_core/opengate_lib/GateUniformElectroMagneticField.h create mode 100644 core/opengate_core/opengate_lib/pyGateElectroMagneticField.cpp create mode 100644 core/opengate_core/opengate_lib/pyGateUniformElectroMagneticField.cpp diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index f09ff7a458..1ba8a32efa 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -211,6 +211,8 @@ void init_G4MagInt_Driver(py::module &); void init_G4ChordFinder(py::module &); void init_GateMagneticField(py::module &); +void init_GateElectroMagneticField(py::module &); +void init_GateUniformElectroMagneticField(py::module &); // void init_GateMappedMagneticField(py::module &); @@ -589,6 +591,8 @@ PYBIND11_MODULE(opengate_core, m) { init_G4MagInt_Driver(m); init_G4ChordFinder(m); init_GateMagneticField(m); + init_GateElectroMagneticField(m); + init_GateUniformElectroMagneticField(m); // init_GateMappedMagneticField(m); init_G4Box(m); diff --git a/core/opengate_core/opengate_lib/GateElectroMagneticField.cpp b/core/opengate_core/opengate_lib/GateElectroMagneticField.cpp new file mode 100644 index 0000000000..0483ba2596 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateElectroMagneticField.cpp @@ -0,0 +1,42 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include "GateElectroMagneticField.h" + +GateElectroMagneticField::GateElectroMagneticField( + G4ElectroMagneticField *inner, const G4VSolid *solid, + std::vector translations, + std::vector rotations, double deltaChordMM) + : G4ElectroMagneticField(), GateField(solid, std::move(translations), + std::move(rotations), deltaChordMM), + m_inner(inner) {} + +void GateElectroMagneticField::GetFieldValue(const G4double Point[4], + G4double *BEfield) const { + const G4ThreeVector worldPoint(Point[0], Point[1], Point[2]); + + const G4AffineTransform *transform = nullptr; + const G4ThreeVector localPoint = + findContainingPlacement(worldPoint, transform); + const G4double localPos[4] = {localPoint.x(), localPoint.y(), localPoint.z(), + Point[3]}; + + G4double localBE[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + m_inner->GetFieldValue(localPos, localBE); + + const G4ThreeVector worldB = + rotateToWorld({localBE[0], localBE[1], localBE[2]}, *transform); + const G4ThreeVector worldE = + rotateToWorld({localBE[3], localBE[4], localBE[5]}, *transform); + + BEfield[0] = worldB.x(); + BEfield[1] = worldB.y(); + BEfield[2] = worldB.z(); + BEfield[3] = worldE.x(); + BEfield[4] = worldE.y(); + BEfield[5] = worldE.z(); +} diff --git a/core/opengate_core/opengate_lib/GateElectroMagneticField.h b/core/opengate_core/opengate_lib/GateElectroMagneticField.h new file mode 100644 index 0000000000..d48ab7f5e6 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateElectroMagneticField.h @@ -0,0 +1,41 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateElectroMagneticField_h +#define GateElectroMagneticField_h + +#include "G4ElectroMagneticField.hh" +#include "G4RotationMatrix.hh" +#include "G4ThreeVector.hh" +#include "GateField.h" +#include + +class G4VSolid; + +// GATE wrapper for G4ElectroMagneticField +class GateElectroMagneticField : public G4ElectroMagneticField, + protected GateField { +public: + GateElectroMagneticField(G4ElectroMagneticField *inner, const G4VSolid *solid, + std::vector translations, + std::vector rotations, + double deltaChordMM); + + void GetFieldValue(const G4double Point[4], G4double *BEfield) const override; + + G4bool DoesFieldChangeEnergy() const override { return true; } + + inline void SetTransforms(std::vector translations, + std::vector rotations) { + GateField::SetTransforms(std::move(translations), std::move(rotations)); + } + +private: + G4ElectroMagneticField *m_inner; +}; + +#endif // GateElectroMagneticField_h diff --git a/core/opengate_core/opengate_lib/GateUniformElectroMagneticField.cpp b/core/opengate_core/opengate_lib/GateUniformElectroMagneticField.cpp new file mode 100644 index 0000000000..fe5c0aa6f6 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateUniformElectroMagneticField.cpp @@ -0,0 +1,24 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include "GateUniformElectroMagneticField.h" + +GateUniformElectroMagneticField::GateUniformElectroMagneticField( + G4ThreeVector e_field_vector, G4ThreeVector b_field_vector) + : G4ElectroMagneticField(), m_e_field_vector(e_field_vector), + m_b_field_vector(b_field_vector) {} + +void GateUniformElectroMagneticField::GetFieldValue(const G4double Point[4], + G4double *BEfield) const { + // Uniform field does not depend on position or time. + BEfield[0] = m_b_field_vector.x(); + BEfield[1] = m_b_field_vector.y(); + BEfield[2] = m_b_field_vector.z(); + BEfield[3] = m_e_field_vector.x(); + BEfield[4] = m_e_field_vector.y(); + BEfield[5] = m_e_field_vector.z(); +} diff --git a/core/opengate_core/opengate_lib/GateUniformElectroMagneticField.h b/core/opengate_core/opengate_lib/GateUniformElectroMagneticField.h new file mode 100644 index 0000000000..bd24bc106c --- /dev/null +++ b/core/opengate_core/opengate_lib/GateUniformElectroMagneticField.h @@ -0,0 +1,32 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateUniformElectroMagneticField_h +#define GateUniformElectroMagneticField_h + +#include "G4ElectroMagneticField.hh" +#include "G4ThreeVector.hh" + +// Pure inner field class (Geant4-like): uniform electric and magnetic fields. +// This does NOT inherit from GateElectroMagneticField; it is a pure Geant4 +// inner field that will be wrapped by GateElectroMagneticField in the Python +// field manager. +class GateUniformElectroMagneticField : public G4ElectroMagneticField { +public: + GateUniformElectroMagneticField(G4ThreeVector e_field_vector, + G4ThreeVector b_field_vector); + + void GetFieldValue(const G4double Point[4], G4double *BEfield) const override; + + G4bool DoesFieldChangeEnergy() const override { return true; } + +private: + G4ThreeVector m_e_field_vector; + G4ThreeVector m_b_field_vector; +}; + +#endif // GateUniformElectroMagneticField_h \ No newline at end of file diff --git a/core/opengate_core/opengate_lib/pyGateElectroMagneticField.cpp b/core/opengate_core/opengate_lib/pyGateElectroMagneticField.cpp new file mode 100644 index 0000000000..67f919fc0a --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateElectroMagneticField.cpp @@ -0,0 +1,36 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include +#include + +namespace py = pybind11; + +#include "G4ElectroMagneticField.hh" +#include "G4VSolid.hh" +#include "GateElectroMagneticField.h" + +void init_GateElectroMagneticField(py::module &m) { + py::class_>( + m, "GateElectroMagneticField") + + .def(py::init([](G4ElectroMagneticField *inner, const G4VSolid *solid, + std::vector translations, + std::vector rotations, + double delta_chord_mm) { + return new GateElectroMagneticField(inner, solid, translations, + rotations, delta_chord_mm); + }), + py::arg("inner_field"), py::arg("solid"), py::arg("translations"), + py::arg("rotations"), py::arg("delta_chord_mm")) + + .def("SetTransforms", &GateElectroMagneticField::SetTransforms, + py::arg("translations"), py::arg("rotations"), + "Replace the cached world-to-local transforms. Call between runs " + "only after dynamic geometry changes have been applied."); +} diff --git a/core/opengate_core/opengate_lib/pyGateUniformElectroMagneticField.cpp b/core/opengate_core/opengate_lib/pyGateUniformElectroMagneticField.cpp new file mode 100644 index 0000000000..b4e55c4e8a --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateUniformElectroMagneticField.cpp @@ -0,0 +1,27 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include +#include + +namespace py = pybind11; + +#include "G4ElectroMagneticField.hh" +#include "GateUniformElectroMagneticField.h" + +void init_GateUniformElectroMagneticField(py::module &m) { + py::class_>( + m, "GateUniformElectroMagneticField") + + .def(py::init( + [](G4ThreeVector e_field_vector, G4ThreeVector b_field_vector) { + return new GateUniformElectroMagneticField(e_field_vector, + b_field_vector); + }), + py::arg("e_field_vector"), py::arg("b_field_vector")); +} diff --git a/opengate/engines.py b/opengate/engines.py index a5f8b7b0fd..ff3f29d580 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -955,6 +955,7 @@ def ConstructSDandField(self): volume_obj = self.volume_manager.get_volume(volume_name) field._field_volume_obj = volume_obj volume_obj.g4_field_manager = field.create_field_manager() + field._g4_runtime_objects[-1]["volume"] = volume_obj volume_obj.g4_logical_volume.SetFieldManager( volume_obj.g4_field_manager, True ) diff --git a/opengate/geometry/fields.py b/opengate/geometry/fields.py index 0b7a6b67a4..d70956596a 100644 --- a/opengate/geometry/fields.py +++ b/opengate/geometry/fields.py @@ -86,6 +86,12 @@ def __init__(self, *args, **kwargs) -> None: self.attached_to = [] self._field_changes_energy = False + # Integration objects — shared by all field subclasses + self.g4_equation_of_motion = None + self.g4_integrator_stepper = None + self.g4_integration_driver = None + self.g4_chord_finder = None + @property def field_changes_energy(self) -> bool: """Whether the field changes particle energy (False for magnetic, True for others).""" @@ -95,102 +101,75 @@ def refresh_transforms(self) -> None: """Recompute and push cached world-to-local transforms after dynamic geometry changes. """ - if self.g4_field is None or self._field_volume_obj is None: - return - - volume = self._field_volume_obj - g4_translations = [] - g4_rotations = [] - for i in range(volume.number_of_repetitions): - pv = volume.get_g4_physical_volume(i) - T = vec_g4_as_np(pv.GetObjectTranslation()) - R = rot_g4_as_np(pv.GetObjectRotation()) - for anc in volume.ancestor_volumes[::-1]: - anc_pv = getattr(anc, "g4_physical_volume", None) - if anc_pv is None: - continue - anc_T = vec_g4_as_np(anc_pv.GetObjectTranslation()) - anc_R = rot_g4_as_np(anc_pv.GetObjectRotation()) - T = anc_R @ T + anc_T - R = anc_R @ R - g4_translations.append(vec_np_as_g4(T)) - g4_rotations.append(rot_np_as_g4(R)) - - self.g4_field.SetTransforms(g4_translations, g4_rotations) + for entry in self._g4_runtime_objects: + volume = entry.get("volume") + g4_field = entry.get("field") + if volume is None or g4_field is None: + continue + + g4_translations = [] + g4_rotations = [] + for i in range(volume.number_of_repetitions): + pv = volume.get_g4_physical_volume(i) + T = vec_g4_as_np(pv.GetObjectTranslation()) + R = rot_g4_as_np(pv.GetObjectRotation()) + for anc in volume.ancestor_volumes[::-1]: + anc_pv = getattr(anc, "g4_physical_volume", None) + if anc_pv is None: + continue + anc_T = vec_g4_as_np(anc_pv.GetObjectTranslation()) + anc_R = rot_g4_as_np(anc_pv.GetObjectRotation()) + T = anc_R @ T + anc_T + R = anc_R @ R + g4_translations.append(vec_np_as_g4(T)) + g4_rotations.append(rot_np_as_g4(R)) + + g4_field.SetTransforms(g4_translations, g4_rotations) def create_field_manager(self) -> g4.G4FieldManager: """Construct the field and return a configured G4FieldManager.""" msg = "create_field_manager() must be implemented in subclasses." raise NotImplementedError(msg) - def to_dictionary(self): - d = super().to_dictionary() - d["attached_to"] = list(self.attached_to) - return d - - def from_dictionary(self, d): - super().from_dictionary(d) - self.attached_to = d.get("attached_to", []) - - def close(self) -> None: - self.g4_field = None - self._g4_runtime_objects = [] - super().close() - - -class MagneticField(FieldBase): - """Base class for magnetic fields.""" - - def __init__(self, *args, **kwargs) -> None: - super().__init__(*args, **kwargs) - self._field_changes_energy = False - self.g4_equation_of_motion = None - self.g4_integrator_stepper = None - self.g4_integration_driver = None - self.g4_chord_finder = None - - def _create_inner_field(self): - """Create and return the inner G4MagneticField in local volume coordinates. - Override in subclasses.""" - raise NotImplementedError( - "_create_inner_field() must be implemented in subclasses." - ) - - def create_field_manager(self) -> g4.G4FieldManager: - """Construct the field and return a configured G4FieldManager.""" - # Build the inner field (subclass responsibility, in local volume coords). - inner = self._create_inner_field() - - # Wrap in GateMagneticField: handles multi-placement coordinate transforms - # and the integration-overshoot fallback uniformly for all field types. + def _make_g4_transforms(self): + """Return (g4_translations, g4_rotations) for the current field volume.""" translations_np, rotations_np = get_transform_world_to_local( self._field_volume_obj ) - self.g4_field = g4.GateMagneticField( - inner, - self._field_volume_obj.g4_solid, + return ( [vec_np_as_g4(t) for t in translations_np], [rot_np_as_g4(r) for r in rotations_np], - self.delta_chord, ) - # TODO: allow user choice of stepper and equation type - self.g4_equation_of_motion = g4.G4Mag_UsualEqRhs(self.g4_field) + @staticmethod + def _validate_field_function(func, class_name, n_components): + """Validate that func is callable and returns exactly n_components values.""" + if func is None: + raise ValueError(f"field_function must be provided for {class_name}") + if not callable(func): + raise TypeError("field_function must be a callable function") + result = list(func(0, 0, 0, 0)) + if len(result) != n_components: + raise ValueError( + f"{class_name}: field_function must return {n_components} components, " + f"got {len(result)}" + ) + + def _build_field_manager(self, inner_field, gate_field, equation_cls, n_vars): + """Build equation/stepper/driver/chord_finder/fm, record runtime objects, return fm.""" + self.g4_field = gate_field + self.g4_equation_of_motion = equation_cls(gate_field) self.g4_integrator_stepper = g4.G4ClassicalRK4( - self.g4_equation_of_motion, - 6, # x,y,z + px,py,pz + self.g4_equation_of_motion, n_vars ) self.g4_integration_driver = g4.G4MagInt_Driver( - self.step_minimum, - self.g4_integrator_stepper, - 6, - 0, # no verbosity + self.step_minimum, self.g4_integrator_stepper, n_vars, 0 ) self.g4_chord_finder = g4.G4ChordFinder(self.g4_integration_driver) self.g4_chord_finder.SetDeltaChord(self.delta_chord) fm = g4.G4FieldManager( - self.g4_field, self.g4_chord_finder, self.field_changes_energy + gate_field, self.g4_chord_finder, self.field_changes_energy ) fm.SetDeltaOneStep(self.delta_one_step) fm.SetDeltaIntersection(self.delta_intersection) @@ -198,12 +177,10 @@ def create_field_manager(self) -> g4.G4FieldManager: fm.SetMaximumEpsilonStep(self.max_epsilon_step) # Keep all runtime objects alive for the full simulation lifetime. - # 'inner_field' is held as a raw pointer by GateMagneticField in C++, - # so Python must retain a reference here. self._g4_runtime_objects.append( { - "inner_field": inner, - "field": self.g4_field, + "inner_field": inner_field, + "field": gate_field, "equation": self.g4_equation_of_motion, "stepper": self.g4_integrator_stepper, "driver": self.g4_integration_driver, @@ -214,14 +191,50 @@ def create_field_manager(self) -> g4.G4FieldManager: return fm + def to_dictionary(self): + d = super().to_dictionary() + d["attached_to"] = list(self.attached_to) + return d + + def from_dictionary(self, d): + super().from_dictionary(d) + self.attached_to = d.get("attached_to", []) + def close(self) -> None: + self.g4_field = None self.g4_chord_finder = None self.g4_integration_driver = None self.g4_integrator_stepper = None self.g4_equation_of_motion = None + self._g4_runtime_objects = [] super().close() +class MagneticField(FieldBase): + """Base class for magnetic fields.""" + + def _create_inner_field(self): + """Create and return the inner G4MagneticField in local volume coordinates. + Override in subclasses.""" + raise NotImplementedError( + "_create_inner_field() must be implemented in subclasses." + ) + + def create_field_manager(self) -> g4.G4FieldManager: + """Construct the field and return a configured G4FieldManager.""" + inner = self._create_inner_field() + g4_translations, g4_rotations = self._make_g4_transforms() + gate_field = g4.GateMagneticField( + inner, + self._field_volume_obj.g4_solid, + g4_translations, + g4_rotations, + self.delta_chord, + ) + # TODO: allow user choice of stepper and equation type + return self._build_field_manager(inner, gate_field, g4.G4Mag_UsualEqRhs, 6) + + class UniformMagneticField(MagneticField): """Uniform magnetic field with constant field vector. @@ -243,9 +256,6 @@ class UniformMagneticField(MagneticField): ), } - def __init__(self, *args, **kwargs) -> None: - super().__init__(*args, **kwargs) - def _create_inner_field(self): return g4.G4UniformMagField(g4.G4ThreeVector(*self.field_vector)) @@ -265,9 +275,6 @@ class QuadrupoleMagneticField(MagneticField): ), } - def __init__(self, *args, **kwargs) -> None: - super().__init__(*args, **kwargs) - def _create_inner_field(self): return g4.G4QuadrupoleMagField(self.gradient) @@ -287,9 +294,6 @@ class CustomMagneticField(MagneticField): ), } - def __init__(self, *args, **kwargs) -> None: - super().__init__(*args, **kwargs) - def _create_inner_field(self): """Create the custom magnetic field using the Python trampoline. @@ -297,17 +301,7 @@ def _create_inner_field(self): the attached volume and must return [Bx, By, Bz] in that same local frame. The base class rotates the result to world coordinates. """ - if self.field_function is None: - raise ValueError("field_function must be provided for CustomMagneticField") - if not callable(self.field_function): - raise TypeError("field_function must be a callable function") - - test_point = [0, 0, 0, 0] - result = list(self.field_function(*test_point)) - if len(result) != 3: - raise ValueError( - "field_function must return a list of 3 components: [Bx, By, Bz]" - ) + self._validate_field_function(self.field_function, "CustomMagneticField", 3) class _PyMagneticField(g4.G4MagneticField): def __init__(inner_self, callback): @@ -326,82 +320,40 @@ def to_dictionary(self): class ElectroMagneticField(FieldBase): - """Base class for electromagnetic fields.""" + """Base class for electromagnetic fields (includes pure electric).""" def __init__(self, *args, **kwargs) -> None: super().__init__(*args, **kwargs) # Electromagnetic fields change particle energy self._field_changes_energy = True - self.g4_equation_of_motion = None - self.g4_integrator_stepper = None - self.g4_integration_driver = None - self.g4_chord_finder = None - - def _create_field(self) -> None: - """Create the G4 field object. Override in subclasses.""" - msg = "_create_field() must be implemented in subclasses." - raise NotImplementedError(msg) + def _create_inner_field(self): + """Create and return the inner G4ElectroMagneticField in local coordinates. + Override in subclasses.""" + raise NotImplementedError( + "_create_inner_field() must be implemented in subclasses." + ) def create_field_manager(self) -> g4.G4FieldManager: """Construct the field and return a configured G4FieldManager.""" - # Create the field (subclass responsibility) - self._create_field() - - self.g4_equation_of_motion = g4.G4EqMagElectricField(self.g4_field) - self.g4_integrator_stepper = g4.G4ClassicalRK4( - self.g4_equation_of_motion, - 8, # number of variables for electromagnetic field = 8 (x,y,z + px,py,pz + t + E) - ) - - self.g4_integration_driver = g4.G4MagInt_Driver( - self.step_minimum, - self.g4_integrator_stepper, - 8, # number of variables for electromagnetic field = 8 (x,y,z + px,py,pz + t + E) - 0, - ) - - self.g4_chord_finder = g4.G4ChordFinder(self.g4_integration_driver) - self.g4_chord_finder.SetDeltaChord(self.delta_chord) - - # Create and configure field manager - fm = g4.G4FieldManager( - self.g4_field, self.g4_chord_finder, self.field_changes_energy - ) - fm.SetDeltaOneStep(self.delta_one_step) - fm.SetDeltaIntersection(self.delta_intersection) - fm.SetMinimumEpsilonStep(self.min_epsilon_step) - fm.SetMaximumEpsilonStep(self.max_epsilon_step) - - # Keep runtime objects alive for the full simulation lifetime. - self._g4_runtime_objects.append( - { - "field": self.g4_field, - "equation": self.g4_equation_of_motion, - "stepper": self.g4_integrator_stepper, - "driver": self.g4_integration_driver, - "chord_finder": self.g4_chord_finder, - "field_manager": fm, - } + inner = self._create_inner_field() + g4_translations, g4_rotations = self._make_g4_transforms() + gate_field = g4.GateElectroMagneticField( + inner, + self._field_volume_obj.g4_solid, + g4_translations, + g4_rotations, + self.delta_chord, ) - - return fm - - def close(self) -> None: - self.g4_chord_finder = None - self.g4_integration_driver = None - self.g4_integrator_stepper = None - self.g4_equation_of_motion = None - super().close() + # TODO: allow user choice of stepper and equation type + # 8 variables: x,y,z + px,py,pz + t + E + return self._build_field_manager(inner, gate_field, g4.G4EqMagElectricField, 8) class ElectricField(ElectroMagneticField): """Base class for pure electric fields.""" - def __init__(self, *args, **kwargs) -> None: - super().__init__(*args, **kwargs) - # Electric fields change particle energy - self._field_changes_energy = True + # _field_changes_energy is already True from ElectroMagneticField class UniformElectricField(ElectricField): @@ -414,17 +366,14 @@ class UniformElectricField(ElectricField): "field_vector": ( [0, 0, 0], { - "doc": "Field vector [Ex, Ey, Ez] in Geant4 units.", + "doc": "Field vector [Ex, Ey, Ez] in local volume coordinates. " + "Each component in electric field strength units.", }, ), } - def __init__(self, *args, **kwargs) -> None: - super().__init__(*args, **kwargs) - - def _create_field(self) -> None: - """Create the uniform electric field.""" - self.g4_field = g4.G4UniformElectricField(g4.G4ThreeVector(*self.field_vector)) + def _create_inner_field(self): + return g4.G4UniformElectricField(g4.G4ThreeVector(*self.field_vector)) class CustomElectricField(ElectricField): @@ -442,25 +391,9 @@ class CustomElectricField(ElectricField): ), } - def __init__(self, *args, **kwargs) -> None: - super().__init__(*args, **kwargs) - - def _create_field(self) -> None: - """Create the custom electric field using the Python trampoline.""" - if self.field_function is None: - raise ValueError("field_function must be provided for CustomElectricField") - if not callable(self.field_function): - raise TypeError("field_function must be a callable function") - - # Check if the function returns 3 components - test_point = [0, 0, 0, 0] - result = self.field_function(*test_point) - if len(result) != 3: - raise ValueError( - "field_function must return a list of 3 components: [Ex, Ey, Ez]" - ) + def _create_inner_field(self): + self._validate_field_function(self.field_function, "CustomElectricField", 3) - # Create a custom G4ElectricField subclass that calls our Python function class _PyElectricField(g4.G4ElectricField): def __init__(inner_self, callback): super().__init__() @@ -472,7 +405,7 @@ def GetFieldValue(inner_self, point): def DoesFieldChangeEnergy(inner_self): return True - self.g4_field = _PyElectricField(self.field_function) + return _PyElectricField(self.field_function) def to_dictionary(self): raise NotImplementedError( @@ -480,7 +413,6 @@ def to_dictionary(self): ) -# TODO (@srmarcballestero): this implementation always goes through the Python trampoline, which compromises performance. Need to implement a more efficient way on the C++ side. class UniformElectroMagneticField(ElectroMagneticField): """Uniform electromagnetic field with constant magnetic and electric field vectors.""" @@ -492,48 +424,22 @@ class UniformElectroMagneticField(ElectroMagneticField): "magnetic_field_vector": ( [0, 0, 0], { - "doc": "Magnetic field vector [Bx, By, Bz] in Geant4 units.", + "doc": "Magnetic field vector [Bx, By, Bz] in local volume coordinates.", }, ), "electric_field_vector": ( [0, 0, 0], { - "doc": "Electric field vector [Ex, Ey, Ez] in Geant4 units.", + "doc": "Electric field vector [Ex, Ey, Ez] in local volume coordinates.", }, ), } - def __init__(self, *args, **kwargs) -> None: - raise NotImplementedError("UniformElectroMagneticField is not implemented yet.") - super().__init__(*args, **kwargs) - - def _create_field(self) -> None: - """Create the uniform electromagnetic field using the Python trampoline.""" - bx, by, bz = self.magnetic_field_vector - ex, ey, ez = self.electric_field_vector - - # Create a custom G4ElectroMagneticField subclass with constant field values - class _PyUniformEMField(g4.G4ElectroMagneticField): - def __init__(inner_self, b_field, e_field): - super().__init__() - inner_self._b_field = b_field - inner_self._e_field = e_field - - def GetFieldValue(inner_self, point): - # Return constant [Bx, By, Bz, Ex, Ey, Ez] - return [ - inner_self._b_field[0], - inner_self._b_field[1], - inner_self._b_field[2], - inner_self._e_field[0], - inner_self._e_field[1], - inner_self._e_field[2], - ] - - def DoesFieldChangeEnergy(inner_self): - return True - - self.g4_field = _PyUniformEMField([bx, by, bz], [ex, ey, ez]) + def _create_inner_field(self): + return g4.GateUniformElectroMagneticField( + g4.G4ThreeVector(*self.electric_field_vector), + g4.G4ThreeVector(*self.magnetic_field_vector), + ) class CustomElectroMagneticField(ElectroMagneticField): @@ -551,27 +457,11 @@ class CustomElectroMagneticField(ElectroMagneticField): ), } - def __init__(self, *args, **kwargs) -> None: - super().__init__(*args, **kwargs) - - def _create_field(self) -> None: - """Create the custom electromagnetic field using the Python trampoline.""" - if self.field_function is None: - raise ValueError( - "field_function must be provided for CustomElectroMagneticField" - ) - if not callable(self.field_function): - raise TypeError("field_function must be a callable function") - - # Check if the function returns 6 components - test_point = [0, 0, 0, 0] - result = self.field_function(*test_point) - if len(result) != 6: - raise ValueError( - "field_function must return a list of 6 components: [Bx, By, Bz, Ex, Ey, Ez]" - ) + def _create_inner_field(self): + self._validate_field_function( + self.field_function, "CustomElectroMagneticField", 6 + ) - # Create a custom G4ElectroMagneticField subclass that calls our Python function class _PyElectroMagneticField(g4.G4ElectroMagneticField): def __init__(inner_self, callback): super().__init__() @@ -583,7 +473,7 @@ def GetFieldValue(inner_self, point): def DoesFieldChangeEnergy(inner_self): return True - self.g4_field = _PyElectroMagneticField(self.field_function) + return _PyElectroMagneticField(self.field_function) def to_dictionary(self): raise NotImplementedError( From 5a1d3e1f7c785f49a11d1a095e0b30718da5b9cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Tue, 5 May 2026 17:49:42 +0200 Subject: [PATCH 05/28] Bind G4SextupoleMagField and add support for it on the Python side --- .../g4_bindings/pyG4SextupoleMagField.cpp | 24 +++++++++++++++++++ core/opengate_core/opengate_core.cpp | 3 +++ opengate/geometry/fields.py | 19 +++++++++++++++ 3 files changed, 46 insertions(+) create mode 100644 core/opengate_core/g4_bindings/pyG4SextupoleMagField.cpp diff --git a/core/opengate_core/g4_bindings/pyG4SextupoleMagField.cpp b/core/opengate_core/g4_bindings/pyG4SextupoleMagField.cpp new file mode 100644 index 0000000000..fe9201cedb --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4SextupoleMagField.cpp @@ -0,0 +1,24 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "G4MagneticField.hh" +#include "G4SextupoleMagField.hh" + +void init_G4SextupoleMagField(py::module &m) { + + py::class_>( + m, "G4SextupoleMagField") + + .def(py::init()) + + ; +} diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index 1ba8a32efa..29b274f948 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -188,6 +188,8 @@ void init_G4UniformMagField(py::module &); void init_G4QuadrupoleMagField(py::module &); +void init_G4SextupoleMagField(py::module &); + void init_G4UniformElectricField(py::module &); void init_G4EquationOfMotion(py::module &); @@ -579,6 +581,7 @@ PYBIND11_MODULE(opengate_core, m) { init_G4ElectricField(m); init_G4UniformMagField(m); init_G4QuadrupoleMagField(m); + init_G4SextupoleMagField(m); init_G4UniformElectricField(m); init_G4EquationOfMotion(m); init_G4Mag_EqRhs(m); diff --git a/opengate/geometry/fields.py b/opengate/geometry/fields.py index d70956596a..f723adadc9 100644 --- a/opengate/geometry/fields.py +++ b/opengate/geometry/fields.py @@ -279,6 +279,25 @@ def _create_inner_field(self): return g4.G4QuadrupoleMagField(self.gradient) +class SextupoleMagneticField(MagneticField): + """Sextupole magnetic field with gradient.""" + + # hints for IDE + gradient: float + + user_info_defaults = { + "gradient": ( + 0, + { + "doc": "Field gradient in magnetic field strength units per unit length (e.g., Tesla/meter).", + }, + ), + } + + def _create_inner_field(self): + return g4.G4SextupoleMagField(self.gradient) + + class CustomMagneticField(MagneticField): """Custom magnetic field defined by a Python callback function.""" From 68f6b4f45d786ea6a9db2b27cfd34e7182dcfa43 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Tue, 5 May 2026 17:49:42 +0200 Subject: [PATCH 06/28] Bind G4SextupoleMagField and add support for it on the Python side --- opengate/geometry/fields.py | 1 - 1 file changed, 1 deletion(-) diff --git a/opengate/geometry/fields.py b/opengate/geometry/fields.py index f723adadc9..3c3b7c622f 100644 --- a/opengate/geometry/fields.py +++ b/opengate/geometry/fields.py @@ -14,7 +14,6 @@ # ! ======= KNOWN TODO'S ======== # ! - implement the possibility of choosing the stepper type and equation type -# ! - bind the sextupole magnetic field geant4 implementation # ! WIP - implement mapped fields (e.g., from a CSV file) # ! - Overhead for custom fields implementation: every GetFieldValue call crosses c++ -> Python -> c++, acquiring the GIL each time, # ! which is very inefficient. Need to implement a more efficient way on the C++ side. From 28642945160348a121fbc111e87483760943f92f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Tue, 5 May 2026 19:36:22 +0200 Subject: [PATCH 07/28] Fixes: - copilot code review. - register SextupoleMagneticField in field_types and process_cls --- opengate/engines.py | 4 +--- opengate/geometry/fields.py | 27 ++++++++++++++++++--------- 2 files changed, 19 insertions(+), 12 deletions(-) diff --git a/opengate/engines.py b/opengate/engines.py index ff3f29d580..a7aff20502 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -953,9 +953,7 @@ def ConstructSDandField(self): for field in self.simulation_engine.simulation.volume_manager.fields.values(): for volume_name in field.attached_to: volume_obj = self.volume_manager.get_volume(volume_name) - field._field_volume_obj = volume_obj - volume_obj.g4_field_manager = field.create_field_manager() - field._g4_runtime_objects[-1]["volume"] = volume_obj + volume_obj.g4_field_manager = field.create_field_manager(volume_obj) volume_obj.g4_logical_volume.SetFieldManager( volume_obj.g4_field_manager, True ) diff --git a/opengate/geometry/fields.py b/opengate/geometry/fields.py index 3c3b7c622f..870d6745b5 100644 --- a/opengate/geometry/fields.py +++ b/opengate/geometry/fields.py @@ -78,9 +78,7 @@ def __init__(self, *args, **kwargs) -> None: self.g4_field = None self._g4_runtime_objects = [] - self._field_volume_obj: Any = ( - None # set by engines.py before create_field_manager() is called - ) + self._field_volume_obj: Any = None self.attached_to = [] self._field_changes_energy = False @@ -125,7 +123,7 @@ def refresh_transforms(self) -> None: g4_field.SetTransforms(g4_translations, g4_rotations) - def create_field_manager(self) -> g4.G4FieldManager: + def create_field_manager(self, volume_obj) -> g4.G4FieldManager: """Construct the field and return a configured G4FieldManager.""" msg = "create_field_manager() must be implemented in subclasses." raise NotImplementedError(msg) @@ -154,7 +152,9 @@ def _validate_field_function(func, class_name, n_components): f"got {len(result)}" ) - def _build_field_manager(self, inner_field, gate_field, equation_cls, n_vars): + def _build_field_manager( + self, inner_field, gate_field, equation_cls, n_vars, volume_obj + ): """Build equation/stepper/driver/chord_finder/fm, record runtime objects, return fm.""" self.g4_field = gate_field self.g4_equation_of_motion = equation_cls(gate_field) @@ -185,6 +185,7 @@ def _build_field_manager(self, inner_field, gate_field, equation_cls, n_vars): "driver": self.g4_integration_driver, "chord_finder": self.g4_chord_finder, "field_manager": fm, + "volume": volume_obj, } ) @@ -219,8 +220,9 @@ def _create_inner_field(self): "_create_inner_field() must be implemented in subclasses." ) - def create_field_manager(self) -> g4.G4FieldManager: + def create_field_manager(self, volume_obj) -> g4.G4FieldManager: """Construct the field and return a configured G4FieldManager.""" + self._field_volume_obj = volume_obj inner = self._create_inner_field() g4_translations, g4_rotations = self._make_g4_transforms() gate_field = g4.GateMagneticField( @@ -231,7 +233,9 @@ def create_field_manager(self) -> g4.G4FieldManager: self.delta_chord, ) # TODO: allow user choice of stepper and equation type - return self._build_field_manager(inner, gate_field, g4.G4Mag_UsualEqRhs, 6) + return self._build_field_manager( + inner, gate_field, g4.G4Mag_UsualEqRhs, 6, volume_obj + ) class UniformMagneticField(MagneticField): @@ -352,8 +356,9 @@ def _create_inner_field(self): "_create_inner_field() must be implemented in subclasses." ) - def create_field_manager(self) -> g4.G4FieldManager: + def create_field_manager(self, volume_obj) -> g4.G4FieldManager: """Construct the field and return a configured G4FieldManager.""" + self._field_volume_obj = volume_obj inner = self._create_inner_field() g4_translations, g4_rotations = self._make_g4_transforms() gate_field = g4.GateElectroMagneticField( @@ -365,7 +370,9 @@ def create_field_manager(self) -> g4.G4FieldManager: ) # TODO: allow user choice of stepper and equation type # 8 variables: x,y,z + px,py,pz + t + E - return self._build_field_manager(inner, gate_field, g4.G4EqMagElectricField, 8) + return self._build_field_manager( + inner, gate_field, g4.G4EqMagElectricField, 8, volume_obj + ) class ElectricField(ElectroMagneticField): @@ -632,6 +639,7 @@ def to_dictionary(self): field_types = { "UniformMagneticField": UniformMagneticField, "QuadrupoleMagneticField": QuadrupoleMagneticField, + "SextupoleMagneticField": SextupoleMagneticField, "CustomMagneticField": CustomMagneticField, # "MappedMagneticField": MappedMagneticField, "UniformElectricField": UniformElectricField, @@ -644,6 +652,7 @@ def to_dictionary(self): process_cls(MagneticField) process_cls(UniformMagneticField) process_cls(QuadrupoleMagneticField) +process_cls(SextupoleMagneticField) process_cls(CustomMagneticField) # process_cls(MappedMagneticField) process_cls(ElectroMagneticField) From 682411dec0bc4f2d2045f6659f3512dcec02fd94 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Wed, 6 May 2026 07:55:45 +0200 Subject: [PATCH 08/28] Fix: robustness of `refresh_transforms` --- opengate/geometry/fields.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/opengate/geometry/fields.py b/opengate/geometry/fields.py index 870d6745b5..abf6d27dfd 100644 --- a/opengate/geometry/fields.py +++ b/opengate/geometry/fields.py @@ -1,5 +1,6 @@ from typing import Any +import numpy as np import opengate_core as g4 from ..base import GateObject, process_cls @@ -109,13 +110,17 @@ def refresh_transforms(self) -> None: for i in range(volume.number_of_repetitions): pv = volume.get_g4_physical_volume(i) T = vec_g4_as_np(pv.GetObjectTranslation()) - R = rot_g4_as_np(pv.GetObjectRotation()) + rot = pv.GetObjectRotation() + R = rot_g4_as_np(rot) if rot is not None else np.eye(3) for anc in volume.ancestor_volumes[::-1]: - anc_pv = getattr(anc, "g4_physical_volume", None) - if anc_pv is None: + # Ancestors are assumed to be singly placed, so we use idx 0 + anc_pvs = getattr(anc, "g4_physical_volumes", None) + if not anc_pvs: continue + anc_pv = anc_pvs[0] anc_T = vec_g4_as_np(anc_pv.GetObjectTranslation()) - anc_R = rot_g4_as_np(anc_pv.GetObjectRotation()) + anc_rot = anc_pv.GetObjectRotation() + anc_R = rot_g4_as_np(anc_rot) if anc_rot is not None else np.eye(3) T = anc_R @ T + anc_T R = anc_R @ R g4_translations.append(vec_np_as_g4(T)) From 6e2e24ca89ef6bba246946d5cb5bdeee95173106 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Wed, 6 May 2026 09:31:49 +0200 Subject: [PATCH 09/28] Enahnce existing tests and add new ones --- .../tests/src/geometry/test099_fields_api.py | 4 +- .../src/geometry/test099_fields_helpers.py | 14 ++ .../test099_fields_multi_volume_refresh.py | 205 ++++++++++++++++++ .../test099_fields_repeated_placements.py | 200 +++++++++++++++++ .../geometry/test099_fields_rotated_volume.py | 169 +++++++++++++++ .../geometry/test099_fields_serialization.py | 18 ++ 6 files changed, 608 insertions(+), 2 deletions(-) create mode 100644 opengate/tests/src/geometry/test099_fields_multi_volume_refresh.py create mode 100644 opengate/tests/src/geometry/test099_fields_repeated_placements.py create mode 100644 opengate/tests/src/geometry/test099_fields_rotated_volume.py diff --git a/opengate/tests/src/geometry/test099_fields_api.py b/opengate/tests/src/geometry/test099_fields_api.py index 294b77e1a6..95bd3f1316 100644 --- a/opengate/tests/src/geometry/test099_fields_api.py +++ b/opengate/tests/src/geometry/test099_fields_api.py @@ -131,7 +131,7 @@ def _expect_error(fn, error_msg): f = fields.CustomMagneticField(name="B_bad") f.field_function = "i'm a string, not a function! will i fool you?" ok = _expect_error( - lambda: f._create_field(), + lambda: f._create_inner_field(), "non-callable field_function", ) is_ok = is_ok and ok @@ -141,7 +141,7 @@ def _expect_error(fn, error_msg): f = fields.CustomMagneticField(name="B_wrong_dims") f.field_function = lambda x, y, z, t: [0, 0, 0, 0, 0, 0] ok = _expect_error( - lambda: f._create_field(), + lambda: f._create_inner_field(), "wrong component count", ) is_ok = is_ok and ok diff --git a/opengate/tests/src/geometry/test099_fields_helpers.py b/opengate/tests/src/geometry/test099_fields_helpers.py index f476178546..f3f75e5ce1 100644 --- a/opengate/tests/src/geometry/test099_fields_helpers.py +++ b/opengate/tests/src/geometry/test099_fields_helpers.py @@ -30,6 +30,20 @@ def cyclotron_radius(T, B, m, q): return p / (q * B * 0.299792458) / g4_m +def magnetic_deflection(T, B, m, q, L): + """ + Expected transverse displacement (chord sagitta) after a particle with + kinetic energy T traverses a uniform magnetic field B over axial depth L. + + Assumes the particle enters parallel to the beam axis and exits through + the back face (i.e. L < cyclotron_radius). + + sagitta = r - sqrt(r^2 - L^2) + """ + r = cyclotron_radius(T, B, m, q) + return r - np.sqrt(r**2 - L**2) + + def build_field_simulation( field_obj, kinetic_energy=10 * gate.g4_units.MeV, diff --git a/opengate/tests/src/geometry/test099_fields_multi_volume_refresh.py b/opengate/tests/src/geometry/test099_fields_multi_volume_refresh.py new file mode 100644 index 0000000000..c11241b913 --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_multi_volume_refresh.py @@ -0,0 +1,205 @@ +#!/usr/bin/env python3 +""" +Test 099 - field attached to two volumes, dynamic geometry between runs. + +Verifies that refresh_transforms() updates the GateMagneticField transform +cache for every volume the field is attached to. + +Setup: + Two boxes share one UniformMagneticField with B along local +Y (3 T). + The boxes are separated in X so each has its own dedicated proton source + and protons do not pass through each other's box. + +box1 is at X = -200 mm and is dynamically rotated 90° about Z in run 1. +box2 is at X = +200 mm and is stationary. + + +Checks: + Run 0: + - box1: deflection in -X, no Y deflection + - box2: deflection in -X, no Y deflection (unchanged by box1 rotation) + + Run 1: + - box1: deflection in -Y, no X deflection + - box2: deflection in -X, no Y deflection (unchanged by box1 rotation) + +""" + +import numpy as np +import uproot +from scipy.spatial.transform import Rotation + +import opengate as gate +from opengate.geometry import fields +from opengate.tests import utility +from test099_fields_helpers import ( + g4_mm, + g4_MeV, + g4_tesla, + g4_eplus, + PROTON_MASS, + magnetic_deflection, +) + +g4_s = gate.g4_units.s + +BOX_SIZE = 100 * g4_mm +SEP = 200 * g4_mm # lateral X separation between box centres +By = 3 * g4_tesla +KE = 10 * g4_MeV + +# Analytical sagitta for a 10 MeV proton in By over BOX_SIZE depth. +EXPECTED_DEFLECTION = magnetic_deflection(KE, By, PROTON_MASS, 1 * g4_eplus, BOX_SIZE) +TOL_DEFLECTED = 0.01 * g4_mm # deflected axis must stay within this of expected +TOL_ZERO = 0.01 * g4_mm # non-deflected axis must stay within this + +VISU = False + + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, output_folder="test099_fields") + + sim = gate.Simulation() + sim.g4_verbose = False + sim.visu = False + sim.random_seed = 99001 + sim.number_of_threads = 1 + sim.output_dir = paths.output + + if VISU: + sim.visu = True + sim.visu_type = "qt" + sim.visu_commands.append("/vis/scene/endOfEventAction accumulate") + sim.visu_commands.append("/vis/scene/add/trajectories smooth") + sim.visu_commands.append("/vis/scene/add/magneticField 20 fullArrow") + + # Two runs: run 0 (box1 identity), run 1 (box1 rotated 90° about Z) + sim.run_timing_intervals = [ + [0, 0.5 * g4_s], + [0.5 * g4_s, 1.0 * g4_s], + ] + + sim.world.size = [1 * gate.g4_units.m, 1 * gate.g4_units.m, 1 * gate.g4_units.m] + sim.world.material = "G4_Galactic" + + # box1: rotates 90° about Z in run 1, placed at X = -SEP + box1 = sim.add_volume("Box", "box1") + box1.size = [BOX_SIZE, BOX_SIZE, BOX_SIZE] + box1.material = "G4_Galactic" + box1.translation = [-SEP, 0, 0] + box1.add_dynamic_parametrisation( + rotation=[ + Rotation.from_euler("z", 0, degrees=True).as_matrix(), # run 0 + Rotation.from_euler("z", 90, degrees=True).as_matrix(), # run 1 + ] + ) + + # box2: stationary control, placed at X = +SEP + box2 = sim.add_volume("Box", "box2") + box2.size = [BOX_SIZE, BOX_SIZE, BOX_SIZE] + box2.material = "G4_Galactic" + box2.translation = [+SEP, 0, 0] + + # Same field on both volumes + field = fields.UniformMagneticField(name="B_uniform") + field.field_vector = [0, By, 0] # B along local +Y + box1.add_field(field) + box2.add_field(field) + + # Source for box1: shoots along +Z, aligned with box1 centre + src1 = sim.add_source("GenericSource", "src_box1") + src1.particle = "proton" + src1.n = [1, 1] + src1.energy.type = "mono" + src1.energy.mono = KE + src1.position.type = "point" + src1.position.translation = [-SEP, 0, -300 * g4_mm] + src1.direction.type = "momentum" + src1.direction.momentum = [0, 0, 1] + + # Source for box2: shoots along +Z, aligned with box2 centre + src2 = sim.add_source("GenericSource", "src_box2") + src2.particle = "proton" + src2.n = [1, 1] + src2.energy.type = "mono" + src2.energy.mono = KE + src2.position.type = "point" + src2.position.translation = [+SEP, 0, -300 * g4_mm] + src2.direction.type = "momentum" + src2.direction.momentum = [0, 0, 1] + + # PhaseSpace on box1 + phsp1 = sim.add_actor("PhaseSpaceActor", "phsp_box1") + phsp1.attached_to = "box1" + phsp1.attributes = ["PostPosition", "RunID"] + phsp1.output_filename = "phsp_box1.root" + phsp1.steps_to_store = "exiting" + phsp1.root_output.write_to_disk = True + + # PhaseSpace on box2 + phsp2 = sim.add_actor("PhaseSpaceActor", "phsp_box2") + phsp2.attached_to = "box2" + phsp2.attributes = ["PostPosition", "RunID"] + phsp2.output_filename = "phsp_box2.root" + phsp2.steps_to_store = "exiting" + phsp2.root_output.write_to_disk = True + + sim.run() + + # --- Read results --- + df1 = uproot.open(str(paths.output / "phsp_box1.root"))["phsp_box1;1"].arrays( + library="pd" + ) + df2 = uproot.open(str(paths.output / "phsp_box2.root"))["phsp_box2;1"].arrays( + library="pd" + ) + + r0b1 = df1[df1["RunID"] == 0].iloc[0] + r1b1 = df1[df1["RunID"] == 1].iloc[0] + r0b2 = df2[df2["RunID"] == 0].iloc[0] + r1b2 = df2[df2["RunID"] == 1].iloc[0] + + # Positions relative to each box's centre (boxes are at X = ±SEP) + x1_r0 = r0b1["PostPosition_X"] - (-SEP) + x1_r1 = r1b1["PostPosition_X"] - (-SEP) + y1_r0 = r0b1["PostPosition_Y"] + y1_r1 = r1b1["PostPosition_Y"] + + x2_r0 = r0b2["PostPosition_X"] - (+SEP) + x2_r1 = r1b2["PostPosition_X"] - (+SEP) + y2_r0 = r0b2["PostPosition_Y"] + y2_r1 = r1b2["PostPosition_Y"] + + print( + f"Analytical expected deflection: {EXPECTED_DEFLECTION:.2f} mm (threshold: {TOL_DEFLECTED:.2f} mm)" + ) + print(f"box1 run0: ΔX={x1_r0:.2f} Y={y1_r0:.2f} mm") + print(f"box1 run1: ΔX={x1_r1:.2f} Y={y1_r1:.2f} mm") + print(f"box2 run0: ΔX={x2_r0:.2f} Y={y2_r0:.2f} mm") + print(f"box2 run1: ΔX={x2_r1:.2f} Y={y2_r1:.2f} mm") + + # Run 0: + ok_r0b1 = ( + np.abs(x1_r0 + EXPECTED_DEFLECTION) < TOL_DEFLECTED and abs(y1_r0) < TOL_ZERO + ) + ok_r0b2 = ( + np.abs(x2_r0 + EXPECTED_DEFLECTION) < TOL_DEFLECTED and abs(y2_r0) < TOL_ZERO + ) + + # Run 1, box2: should be unchanged from run 0 + ok_r1b2 = ( + np.abs(x2_r1 + EXPECTED_DEFLECTION) < TOL_DEFLECTED and abs(y2_r1) < TOL_ZERO + ) + + # Run 1, box1: deflection should now be in -Y, no X deflection + ok_r1b1 = ( + np.abs(y1_r1 + EXPECTED_DEFLECTION) < TOL_DEFLECTED and abs(x1_r1) < TOL_ZERO + ) + + print(f"\nbox1 run0 deflects in -X: {ok_r0b1}") + print(f"box2 run0 deflects in -X: {ok_r0b2}") + print(f"box2 run1 deflects in -X (unchanged): {ok_r1b2}") + print(f"box1 run1 deflects in -Y (verifies refresh): {ok_r1b1}") + + is_ok = ok_r0b1 and ok_r0b2 and ok_r1b2 and ok_r1b1 + utility.test_ok(is_ok) diff --git a/opengate/tests/src/geometry/test099_fields_repeated_placements.py b/opengate/tests/src/geometry/test099_fields_repeated_placements.py new file mode 100644 index 0000000000..29bde8d6f3 --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_repeated_placements.py @@ -0,0 +1,200 @@ +#!/usr/bin/env python3 +""" +Test 099 - field on a repeated volume. + +Two scenarios run in a single simulation. + + LEFT (X = -100 mm): single box, depth BOX_SIZE along Z. + RIGHT (X = +100 mm): N_REPS contiguous slabs, each of depth BOX_SIZE/N_REPS, + total depth = BOX_SIZE. + +Both scenarios expose the proton to the same total path length through field +so the X-deflection should match the analytical sagitta for BOX_SIZE in all cases. + +Checks: + 1. Both deflections match the analytical sagitta within TOL_DEFLECTED. + 2. No Y deflection on either side. + 3. Energy conservation (magnetic field does no work). +""" + +import uproot + +import opengate as gate +from opengate.geometry import fields +from opengate.geometry.utility import get_grid_repetition +from opengate.tests import utility + +from test099_fields_helpers import ( + g4_mm, + g4_MeV, + g4_tesla, + g4_eplus, + PROTON_MASS, + magnetic_deflection, +) + +BOX_SIZE = 50 * g4_mm +N_REPS = 10 +SLAB_SIZE = BOX_SIZE / N_REPS +SPACING = SLAB_SIZE +By = 3 * g4_tesla +KE = 10 * g4_MeV +SEP = 100 * g4_mm + +# Analytical sagitta for a proton traversing BOX_SIZE of field. +EXPECTED_DEFLECTION = magnetic_deflection(KE, By, PROTON_MASS, 1 * g4_eplus, BOX_SIZE) + +TOL_DEFLECTED = 0.1 * g4_mm # comparison tolerance +TOL_Y = 0.1 * g4_mm +TOL_E = 0.01 * g4_MeV + +VISU = False + + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, output_folder="test099_fields") + + sim = gate.Simulation() + sim.g4_verbose = False + sim.visu = False + sim.random_seed = 12345 + sim.number_of_threads = 1 + sim.output_dir = paths.output + + if VISU: + sim.visu = True + sim.visu_type = "qt" + sim.visu_commands.append("/vis/scene/endOfEventAction accumulate") + sim.visu_commands.append("/vis/scene/add/trajectories smooth") + sim.visu_commands.append("/vis/scene/add/magneticField 50 fullArrow") + + m = gate.g4_units.m + sim.world.size = [2 * m, 1 * m, 1 * m] + sim.world.material = "G4_Galactic" + + # --- LEFT: single box at X = -SEP --- + box_single = sim.add_volume("Box", "box_single") + box_single.size = [BOX_SIZE, BOX_SIZE, BOX_SIZE] + box_single.material = "G4_Galactic" + box_single.translation = [-SEP, 0, 0] + + field_single = fields.UniformMagneticField(name="B_single") + field_single.field_vector = [0, By, 0] + box_single.add_field(field_single) + + src_single = sim.add_source("GenericSource", "src_single") + src_single.particle = "proton" + src_single.n = 1 + src_single.energy.type = "mono" + src_single.energy.mono = KE + src_single.position.type = "point" + src_single.position.translation = [-SEP, 0, -200 * g4_mm] + src_single.direction.type = "momentum" + src_single.direction.momentum = [0, 0, 1] + + phsp_single = sim.add_actor("PhaseSpaceActor", "phsp_single") + phsp_single.attached_to = "box_single" + phsp_single.attributes = ["PostKineticEnergy", "PostPosition"] + phsp_single.output_filename = "phsp_single.root" + phsp_single.steps_to_store = "exiting" + phsp_single.root_output.write_to_disk = True + + # RIGHT: repeated box (N_REPS placements along Z) at X = +SEP + translations = get_grid_repetition( + [1, 1, N_REPS], [0, 0, SPACING], start=[0, 0, -SPACING] + ) + translations = [[+SEP + t[0], t[1], t[2]] for t in translations] + + box_multi = sim.add_volume("Box", "box_multi") + box_multi.size = [BOX_SIZE, BOX_SIZE, SLAB_SIZE] + box_multi.material = "G4_Galactic" + box_multi.translation = translations + + field_multi = fields.UniformMagneticField(name="B_multi") + field_multi.field_vector = [0, By, 0] + box_multi.add_field(field_multi) + + src_multi = sim.add_source("GenericSource", "src_multi") + src_multi.particle = "proton" + src_multi.n = 1 + src_multi.energy.type = "mono" + src_multi.energy.mono = KE + src_multi.position.type = "point" + src_multi.position.translation = [+SEP, 0, -200 * g4_mm] + src_multi.direction.type = "momentum" + src_multi.direction.momentum = [0, 0, 1] + + phsp_multi = sim.add_actor("PhaseSpaceActor", "phsp_multi") + phsp_multi.attached_to = "box_multi" + phsp_multi.attributes = ["PostKineticEnergy", "PostPosition"] + phsp_multi.output_filename = "phsp_multi.root" + phsp_multi.steps_to_store = "exiting" + phsp_multi.root_output.write_to_disk = True + + sim.run() + + # Read results + df_single = uproot.open(str(paths.output / "phsp_single.root"))[ + "phsp_single;1" + ].arrays(library="pd") + df_multi = uproot.open(str(paths.output / "phsp_multi.root"))[ + "phsp_multi;1" + ].arrays(library="pd") + + row_single = df_single.sort_values("PostPosition_Z").iloc[-1] + row_multi = df_multi.sort_values("PostPosition_Z").iloc[-1] + + # X position relative to each scenario's X centre + x_single = row_single["PostPosition_X"] - (-SEP) + y_single = row_single["PostPosition_Y"] + ke_single = row_single["PostKineticEnergy"] + + x_multi = row_multi["PostPosition_X"] - (+SEP) + y_multi = row_multi["PostPosition_Y"] + ke_multi = row_multi["PostKineticEnergy"] + + print( + f"Analytical expected deflection: {-EXPECTED_DEFLECTION:.2f} mm (tolerance: ±{TOL_DEFLECTED:.2f} mm)" + ) + print( + f"Single box ({BOX_SIZE:.0f} mm): ΔX={x_single:.2f} mm Y={y_single:.2f} mm KE={ke_single:.4f} MeV" + ) + print( + f"{N_REPS}x slabs ({SLAB_SIZE:.1f} mm each): ΔX={x_multi:.2f} mm Y={y_multi:.2f} mm KE={ke_multi:.4f} MeV" + ) + + # 1. Both deflections match the analytical sagitta. + ok_single = abs(x_single + EXPECTED_DEFLECTION) < TOL_DEFLECTED + ok_multi = abs(x_multi + EXPECTED_DEFLECTION) < TOL_DEFLECTED + print( + f"Single box matches analytical: {ok_single} (|{x_single:.2f} - {-EXPECTED_DEFLECTION:.2f}| < {TOL_DEFLECTED:.2f})" + ) + print( + f"{N_REPS}-slabs match analytical: {ok_multi} (|{x_multi:.2f} - {-EXPECTED_DEFLECTION:.2f}| < {TOL_DEFLECTED:.2f})" + ) + + # 2. No Y deflection + ok_y_single = abs(y_single) < TOL_Y + ok_y_multi = abs(y_multi) < TOL_Y + print(f"No Y deflection (single): {ok_y_single} ({y_single:.3f} mm)") + print(f"No Y deflection (multi): {ok_y_multi} ({y_multi:.3f} mm)") + + # 3. Energy conservation + ok_energy_single = abs(ke_single - KE) < TOL_E + ok_energy_multi = abs(ke_multi - KE) < TOL_E + print( + f"Energy conserved (single): {ok_energy_single} (dKE={ke_single - KE:.4f} MeV)" + ) + print( + f"Energy conserved (multi): {ok_energy_multi} (dKE={ke_multi - KE:.4f} MeV)" + ) + + is_ok = ( + ok_single + and ok_multi + and ok_y_single + and ok_y_multi + and ok_energy_single + and ok_energy_multi + ) + utility.test_ok(is_ok) diff --git a/opengate/tests/src/geometry/test099_fields_rotated_volume.py b/opengate/tests/src/geometry/test099_fields_rotated_volume.py new file mode 100644 index 0000000000..36fe8bfaa2 --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_rotated_volume.py @@ -0,0 +1,169 @@ +#!/usr/bin/env python3 +""" +Test 099 - field direction follows volume rotation. + +Two boxes share one UniformMagneticField with B along local +Y (3 T). +box_a is unrotated, box_b is rotated 90° about Z. + +Checks: + - box_a: deflection in -X, no Y deflection + - box_b: deflection in -Y, no X deflection + - box_a deflection in -X is similar to box_b deflection in -Y +""" + +import numpy as np +import uproot +from scipy.spatial.transform import Rotation + +import opengate as gate +from opengate.geometry import fields +from opengate.tests import utility + +from test099_fields_helpers import ( + g4_mm, + g4_MeV, + g4_tesla, + g4_eplus, + PROTON_MASS, + magnetic_deflection, +) + +BOX_SIZE = 100 * g4_mm +By = 3 * g4_tesla +KE = 10 * g4_MeV +SEP = 300 * g4_mm + +# Expected sagitta for a 10 MeV proton in By over BOX_SIZE depth. +EXPECTED_DEFLECTION = magnetic_deflection(KE, By, PROTON_MASS, 1 * g4_eplus, BOX_SIZE) +TOL_DEFLECTED = 0.01 * g4_mm # deflected axis must stay within this of expected +TOL_ZERO = 0.01 * g4_mm # non-deflected axis must stay within this + +VISU = False # set True for interactive qt viewer + + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, output_folder="test099_fields") + + sim = gate.Simulation() + sim.g4_verbose = False + sim.visu = False + sim.random_seed = 42 + sim.number_of_threads = 1 + sim.output_dir = paths.output + + if VISU: + sim.visu = True + sim.visu_type = "qt" + sim.visu_commands.append("/vis/scene/endOfEventAction accumulate") + sim.visu_commands.append("/vis/scene/add/trajectories smooth") + sim.visu_commands.append("/vis/scene/add/magneticField 20 fullArrow") + + m = gate.g4_units.m + sim.world.size = [2 * m, 2 * m, 2 * m] + sim.world.material = "G4_Galactic" + + # box_a: no rotation + box_a = sim.add_volume("Box", "box_a") + box_a.size = [BOX_SIZE, BOX_SIZE, BOX_SIZE] + box_a.material = "G4_Galactic" + box_a.translation = [-SEP, 0, 0] + + # box_b: rotated 90° around Z + box_b = sim.add_volume("Box", "box_b") + box_b.size = [BOX_SIZE, BOX_SIZE, BOX_SIZE] + box_b.material = "G4_Galactic" + box_b.translation = [+SEP, 0, 0] + box_b.rotation = Rotation.from_euler("z", 90, degrees=True).as_matrix() + + # Same field on both boxes: B along local +Y + field = fields.UniformMagneticField(name="B_shared") + field.field_vector = [0, By, 0] + box_a.add_field(field) + box_b.add_field(field) + + # src_a: shoots through box_a + src_a = sim.add_source("GenericSource", "src_a") + src_a.particle = "proton" + src_a.n = 1 + src_a.energy.type = "mono" + src_a.energy.mono = KE + src_a.position.type = "point" + src_a.position.translation = [-SEP, 0, -300 * g4_mm] + src_a.direction.type = "momentum" + src_a.direction.momentum = [0, 0, 1] + + # src_b: shoots through box_b + src_b = sim.add_source("GenericSource", "src_b") + src_b.particle = "proton" + src_b.n = 1 + src_b.energy.type = "mono" + src_b.energy.mono = KE + src_b.position.type = "point" + src_b.position.translation = [+SEP, 0, -300 * g4_mm] + src_b.direction.type = "momentum" + src_b.direction.momentum = [0, 0, 1] + + # phsp actors + phsp_a = sim.add_actor("PhaseSpaceActor", "phsp_a") + phsp_a.attached_to = "box_a" + phsp_a.attributes = ["PostPosition"] + phsp_a.output_filename = "phsp_a.root" + phsp_a.steps_to_store = "exiting" + phsp_a.root_output.write_to_disk = True + + phsp_b = sim.add_actor("PhaseSpaceActor", "phsp_b") + phsp_b.attached_to = "box_b" + phsp_b.attributes = ["PostPosition"] + phsp_b.output_filename = "phsp_b.root" + phsp_b.steps_to_store = "exiting" + phsp_b.root_output.write_to_disk = True + + sim.run() + + df_a = uproot.open(str(paths.output / "phsp_a.root"))["phsp_a;1"].arrays( + library="pd" + ) + df_b = uproot.open(str(paths.output / "phsp_b.root"))["phsp_b;1"].arrays( + library="pd" + ) + + # Last exit position + row_a = df_a.sort_values("PostPosition_Z").iloc[-1] + row_b = df_b.sort_values("PostPosition_Z").iloc[-1] + + # Positions in the frame of each box centre + xa = row_a["PostPosition_X"] - (-SEP) + ya = row_a["PostPosition_Y"] + xb = row_b["PostPosition_X"] - (+SEP) + yb = row_b["PostPosition_Y"] + + print( + f"Analytical expected deflection: {-EXPECTED_DEFLECTION:.2f} mm (threshold: {TOL_DEFLECTED:.2f} mm)" + ) + print(f"box_a (unrotated): X={xa:.2f} mm Y={ya:.2f} mm") + print(f"box_b (R_z 90°): X={xb:.2f} mm Y={yb:.2f} mm") + + # box_a: deflected in -X, no Y deflection + ok_a_x = np.abs(xa + EXPECTED_DEFLECTION) < TOL_DEFLECTED + ok_a_y = np.abs(ya) < TOL_ZERO + print( + f"box_a deflects in -X: {ok_a_x}, (|{xa:.2f} - {-EXPECTED_DEFLECTION:.2f}| < {TOL_DEFLECTED:.2f})" + ) + print(f"box_a no Y deflection: {ok_a_y} ({ya:.2f} mm)") + + # box_b: B_world = -X → deflects in -Y, no X deflection + ok_b_y = np.abs(yb + EXPECTED_DEFLECTION) < TOL_DEFLECTED + ok_b_x = np.abs(xb) < TOL_ZERO + print( + f"box_b deflects in -Y: {ok_b_y}, (|{yb:.2f} - {-EXPECTED_DEFLECTION:.2f}| < {TOL_DEFLECTED:.2f})" + ) + print(f"box_b no X deflection: {ok_b_x} ({xb:.2f} mm)") + + # box_a deflection in -X is similar to box_b deflection in -Y + ok_a_x_similar_b_y = np.abs(xa - yb) < TOL_DEFLECTED + print( + f"box_a deflection in -X is similar to box_b deflection in -Y: {ok_a_x_similar_b_y}" + ) + + is_ok = ok_a_x and ok_a_y and ok_b_y and ok_b_x and ok_a_x_similar_b_y + utility.test_ok(is_ok) diff --git a/opengate/tests/src/geometry/test099_fields_serialization.py b/opengate/tests/src/geometry/test099_fields_serialization.py index 727bfc09fa..d3c80adb5b 100644 --- a/opengate/tests/src/geometry/test099_fields_serialization.py +++ b/opengate/tests/src/geometry/test099_fields_serialization.py @@ -76,6 +76,24 @@ def _make_sim_with_box(): is_ok = is_ok and ok print(f"Quadrupole round-trip: {'OK' if ok else 'FAIL'}") + # Sextupole round-trip + sim, box = _make_sim_with_box() + field = fields.SextupoleMagneticField(name="B_sext") + field.gradient = 5 * g4_tesla / g4_m + box.add_field(field) + + d = sim.to_dictionary() + sim2 = gate.Simulation() + sim2.from_dictionary(d) + restored = sim2.volume_manager.fields["B_sext"] + ok = ( + isinstance(restored, fields.SextupoleMagneticField) + and restored.gradient == 5 * g4_tesla / g4_m + and restored.attached_to == ["box"] + ) + is_ok = is_ok and ok + print(f"Sextupole round-trip: {'OK' if ok else 'FAIL'}") + # Uniform electric field round-trip sim, box = _make_sim_with_box() field = fields.UniformElectricField(name="E_uniform") From 98d66894a0caea69e6b89864b2334c2de672a232 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Wed, 6 May 2026 09:57:39 +0200 Subject: [PATCH 10/28] Document changes so far --- docs/source/user_guide/user_guide_fields.rst | 56 ++++++++++++++++---- opengate/geometry/fields.py | 6 +-- 2 files changed, 48 insertions(+), 14 deletions(-) diff --git a/docs/source/user_guide/user_guide_fields.rst b/docs/source/user_guide/user_guide_fields.rst index 247bd0d06d..b8cfa8e2d6 100644 --- a/docs/source/user_guide/user_guide_fields.rst +++ b/docs/source/user_guide/user_guide_fields.rst @@ -3,7 +3,11 @@ How to: Electromagnetic fields GATE provides an interface to define electromagnetic fields within simulation volumes. Under the hood, GATE configures Geant4's field tracking machinery (equation of motion, Runge-Kutta stepper, chord finder) for you. -Fields are defined as standalone objects, then attached to one or more volumes with :meth:`~opengate.geometry.volumes.VolumeBase.add_field`. Each volume can hold at most one field. A single field object can be shared across multiple volumes. +Fields are defined as standalone objects, then attached to one or more volumes with ``add_field()``. Each volume can hold at most one field. A single field object can be shared across multiple volumes. + +Fields are attached to logical volumes, so they automatically propagate to all physical placements of that volume. Dynamic geometry changes are supported: if a field is attached to a volume that is later moved or rotated, the field will move/rotate with it. + +All field definitions are relative to the local coordinate system they are attached to. This means that, for instance, if a field is attached to a volume that is rotated, the field vector will be rotated accordingly. Quick example ------------- @@ -37,10 +41,10 @@ Magnetic fields .. code-block:: python field = fields.UniformMagneticField(name="B_uniform") - field.field_vector = [0, 0, 2 * tesla] # [Bx, By, Bz] + field.field_vector = [0, 0, 2 * tesla] # [Bx, By, Bz] <- relative to the volume's local coordinate system box.add_field(field) -**QuadrupoleMagneticField** -- Quadrupole magnetic field defined via a field gradient. Uses the native Geant4 ``G4QuadrupoleMagField`` under the hood. +**QuadrupoleMagneticField** -- Quadrupole magnetic field defined via a field gradient. Uses the native Geant4 ``G4QuadrupoleMagField`` under the hood. The field is always oriented such that the optical axis is along the local Z direction of the volume and one of the south poles sits in the +XY quadrant. .. code-block:: python @@ -48,13 +52,21 @@ Magnetic fields field.gradient = 10 * tesla / m # field gradient (T/m) box.add_field(field) +**SextupoleMagneticField** -- Sextupole magnetic field defined via a field gradient. Uses the native Geant4 ``G4SextupoleMagField`` under the hood. The field is always oriented such that the optical axis is along the local Z direction of the volume and one of the north poles sits in the +Y direction. + +.. code-block:: python + + field = fields.SextupoleMagneticField(name="B_sext") + field.gradient = 10 * tesla / m # field gradient (T/m) + box.add_field(field) + **CustomMagneticField** -- Arbitrary magnetic field defined by a Python callback. The function receives ``(x, y, z, t)`` in Geant4 internal units and must return ``[Bx, By, Bz]``. .. code-block:: python - def my_B_field(x, y, z, t): + def my_B_field(x, y, z, t): # <- relative to the volume's local coordinate system # Spatially varying field - return [0, (1 + x * z / m**2) * tesla, 0] + return [0, (1 + x * z / m**2) * tesla, 0] # <- relative to the volume's local coordinate system field = fields.CustomMagneticField(name="B_custom") field.field_function = my_B_field @@ -76,15 +88,15 @@ Electric fields m = gate.g4_units.m field = fields.UniformElectricField(name="E_uniform") - field.field_vector = [1e6 * volt / m, 0, 0] # [Ex, Ey, Ez] + field.field_vector = [1e6 * volt / m, 0, 0] # [Ex, Ey, Ez] <- relative to the volume's local coordinate system box.add_field(field) **CustomElectricField** -- Arbitrary electric field defined by a Python callback. The function receives ``(x, y, z, t)`` in Geant4 internal units and must return ``[Ex, Ey, Ez]``. .. code-block:: python - def my_E_field(x, y, z, t): - return [1e6 * volt / m, 0, 0] + def my_E_field(x, y, z, t): # <- relative to the volume's local coordinate system + return [1e6 * volt / m, 0, 0] # <- relative to the volume's local coordinate system field = fields.CustomElectricField(name="E_custom") field.field_function = my_E_field @@ -96,6 +108,17 @@ Electromagnetic fields Combined magnetic and electric fields. +**UniformElectroMagneticField** -- Constant electromagnetic field. Uses the GATE C++ implementation ``GateUniformElectroMagneticField`` under the hood. + +.. code-block:: python + + volt = gate.g4_units.volt + m = gate.g4_units.m + + field = fields.UniformElectroMagneticField(name="EM_uniform") + field.field_vector = [0, 0, 1 * tesla, 1e6 * volt / m, 0, 0] # [Bx, By, Bz, Ex, Ey, Ez] <- relative to the volume's local coordinate system + box.add_field(field) + **CustomElectroMagneticField** -- Arbitrary combined field. The callback must return all six components ``[Bx, By, Bz, Ex, Ey, Ez]``. .. code-block:: python @@ -153,7 +176,7 @@ Example: Attaching fields to volumes ---------------------------- -Use :meth:`~opengate.geometry.volumes.VolumeBase.add_field` to attach a field to a volume. The volume must already be added to the simulation. +Use ``add_field()`` to attach a field to a volume. The volume must already be added to the simulation. .. code-block:: python @@ -177,6 +200,10 @@ Use :meth:`~opengate.geometry.volumes.VolumeBase.add_field` to attach a field to **Propagation to daughter volumes.** A field attached to a volume automatically propagates to all of its daughter volumes. A daughter volume can override this by attaching its own field. In that case, the parent's field stops at the daughter's boundary and the daughter's field is used inside. This mirrors standard Geant4 behaviour. +**Behaviour with repeated placements.** When a volume with a field is repeatedly placed (i.e. a logical volume with several physical placements), each physical instance will have the field in its own local coordinate system. + +**Behaviour with dynamic geometry changes.** If a field is attached to a volume that is later moved or rotated, the field will move/rotate with it. + .. code-block:: python outer = sim.add_volume("Box", "outer") @@ -204,7 +231,8 @@ Geant4 can overlay field arrows on the visualization of the geometry. The Geant4 sim.visu = True sim.visu_type = "qt" - sim.visu_commands.append("/vis/scene/add/magneticField 20 lightArrow") + sim.visu_commands.append("/vis/scene/add/magneticField 20 fullArrow") + sim.visu_commands.append("/vis/scene/add/electricField 20 fullArrow") .. _user_guide_fields_performance: @@ -212,7 +240,7 @@ Geant4 can overlay field arrows on the visualization of the geometry. The Geant4 Performance: native vs custom fields ------------------------------------ -Native field types (``UniformMagneticField``, ``UniformElectricField``, ``QuadrupoleMagneticField``) are evaluated entirely in C++ by Geant4. They have no Python overhead and are the recommended choice when they match your use case. +Native field types (``UniformMagneticField``, ``UniformElectricField``, ``QuadrupoleMagneticField``, ``SextupoleMagneticField``, ``UniformElectroMagneticField``) are evaluated entirely in C++ by Geant4 (or GATE in the case of the uniform electromagnetic field). They have no Python overhead and are the recommended choice when they match your use case. Custom fields (``CustomMagneticField``, ``CustomElectricField``, ``CustomElectroMagneticField``) call a Python function for **every evaluation** of ``GetFieldValue`` during tracking. This means that the GIL is acquired and released on every call, which can significantly slow down the simulation, especially in multithreaded mode where all threads serialize through the Python callback. @@ -243,6 +271,9 @@ The field implementation is covered by the ``test099_fields_*`` tests in ``openg - ``test099_fields_analytical_E`` -- Uniform E field vs analytical energy gain. - ``test099_fields_custom_vs_native_B`` -- Custom trampoline B vs native G4 (bit-identical). - ``test099_fields_custom_vs_native_E`` -- Custom trampoline E vs native G4 (bit-identical). +- ``test099_fields_multi_volume_refresh`` -- Testing field updates when volumes are dynamically modified between runs. +- ``test099_fields_repeated_placements`` -- Testing field behaviour with repeated physical placements of the same logical volume. +- ``test099_fields_rotated_volume`` -- Test that fields rotate with their attached volume. - ``test099_fields_serialization`` -- Round-trip serialization for all non-custom types. - ``test099_fields_api`` -- API guards. @@ -259,6 +290,9 @@ Class reference .. autoclass:: opengate.geometry.fields.QuadrupoleMagneticField :members: +.. autoclass:: opengate.geometry.fields.SextupoleMagneticField + :members: + .. autoclass:: opengate.geometry.fields.CustomMagneticField :members: diff --git a/opengate/geometry/fields.py b/opengate/geometry/fields.py index abf6d27dfd..ddebb62188 100644 --- a/opengate/geometry/fields.py +++ b/opengate/geometry/fields.py @@ -316,7 +316,7 @@ class CustomMagneticField(MagneticField): "field_function": ( None, { - "doc": "Python function that takes [x, y, z, t] and returns [Bx, By, Bz].", + "doc": "Python function that takes [x, y, z, t] and returns [Bx, By, Bz], all in local volume coordinates.", }, ), } @@ -416,7 +416,7 @@ class CustomElectricField(ElectricField): "field_function": ( None, { - "doc": "Python function that takes [x, y, z, t] and returns [Ex, Ey, Ez].", + "doc": "Python function that takes [x, y, z, t] and returns [Ex, Ey, Ez], all in local volume coordinates.", }, ), } @@ -482,7 +482,7 @@ class CustomElectroMagneticField(ElectroMagneticField): "field_function": ( None, { - "doc": "Python function that takes [x, y, z, t] and returns [Bx, By, Bz, Ex, Ey, Ez].", + "doc": "Python function that takes [x, y, z, t] and returns [Bx, By, Bz, Ex, Ey, Ez], all in local volume coordinates.", }, ), } From 3f55744a971299aac3ff83f06a446d8eb192ca62 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Wed, 6 May 2026 14:32:36 +0200 Subject: [PATCH 11/28] Refactor internal GATE field management: one single base class contains all the common functionality and field-specific implementations are now minimal Co-authored-by: Copilot --- .../opengate_lib/GateElectroMagneticField.cpp | 32 +++----- .../opengate_lib/GateElectroMagneticField.h | 17 ++-- core/opengate_core/opengate_lib/GateField.h | 61 -------------- .../{GateField.cpp => GateFieldBase.cpp} | 68 +++++++--------- .../opengate_lib/GateFieldBase.h | 81 +++++++++++++++++++ .../opengate_lib/GateMagneticField.cpp | 35 ++------ .../opengate_lib/GateMagneticField.h | 38 +++------ 7 files changed, 142 insertions(+), 190 deletions(-) delete mode 100644 core/opengate_core/opengate_lib/GateField.h rename core/opengate_core/opengate_lib/{GateField.cpp => GateFieldBase.cpp} (55%) create mode 100644 core/opengate_core/opengate_lib/GateFieldBase.h diff --git a/core/opengate_core/opengate_lib/GateElectroMagneticField.cpp b/core/opengate_core/opengate_lib/GateElectroMagneticField.cpp index 0483ba2596..fc143ca6ad 100644 --- a/core/opengate_core/opengate_lib/GateElectroMagneticField.cpp +++ b/core/opengate_core/opengate_lib/GateElectroMagneticField.cpp @@ -11,32 +11,18 @@ GateElectroMagneticField::GateElectroMagneticField( G4ElectroMagneticField *inner, const G4VSolid *solid, std::vector translations, std::vector rotations, double deltaChordMM) - : G4ElectroMagneticField(), GateField(solid, std::move(translations), - std::move(rotations), deltaChordMM), + : G4ElectroMagneticField(), + GateFieldBase(solid, std::move(translations), std::move(rotations), + deltaChordMM), m_inner(inner) {} void GateElectroMagneticField::GetFieldValue(const G4double Point[4], G4double *BEfield) const { - const G4ThreeVector worldPoint(Point[0], Point[1], Point[2]); + // localFieldFunc is a lambda that captures 'this' and calls + // m_inner->GetFieldValue + auto localFieldFunc = [this](const G4double pos[4], G4double *f) { + m_inner->GetFieldValue(pos, f); + }; - const G4AffineTransform *transform = nullptr; - const G4ThreeVector localPoint = - findContainingPlacement(worldPoint, transform); - const G4double localPos[4] = {localPoint.x(), localPoint.y(), localPoint.z(), - Point[3]}; - - G4double localBE[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - m_inner->GetFieldValue(localPos, localBE); - - const G4ThreeVector worldB = - rotateToWorld({localBE[0], localBE[1], localBE[2]}, *transform); - const G4ThreeVector worldE = - rotateToWorld({localBE[3], localBE[4], localBE[5]}, *transform); - - BEfield[0] = worldB.x(); - BEfield[1] = worldB.y(); - BEfield[2] = worldB.z(); - BEfield[3] = worldE.x(); - BEfield[4] = worldE.y(); - BEfield[5] = worldE.z(); + applyTransforms(Point, BEfield, 6, localFieldFunc); } diff --git a/core/opengate_core/opengate_lib/GateElectroMagneticField.h b/core/opengate_core/opengate_lib/GateElectroMagneticField.h index d48ab7f5e6..10a3e22c0a 100644 --- a/core/opengate_core/opengate_lib/GateElectroMagneticField.h +++ b/core/opengate_core/opengate_lib/GateElectroMagneticField.h @@ -9,32 +9,29 @@ #define GateElectroMagneticField_h #include "G4ElectroMagneticField.hh" -#include "G4RotationMatrix.hh" -#include "G4ThreeVector.hh" -#include "GateField.h" +#include "GateFieldBase.h" #include class G4VSolid; -// GATE wrapper for G4ElectroMagneticField +// GATE wrapper for G4ElectroMagneticField. class GateElectroMagneticField : public G4ElectroMagneticField, - protected GateField { + public GateFieldBase { public: + // constructor GateElectroMagneticField(G4ElectroMagneticField *inner, const G4VSolid *solid, std::vector translations, std::vector rotations, double deltaChordMM); + // override GetFieldValue to apply the coordinate transforms void GetFieldValue(const G4double Point[4], G4double *BEfield) const override; + // override DoesFieldChangeEnergy to return true for electro-magnetic fields G4bool DoesFieldChangeEnergy() const override { return true; } - inline void SetTransforms(std::vector translations, - std::vector rotations) { - GateField::SetTransforms(std::move(translations), std::move(rotations)); - } - private: + // the inner field in the volume's local frame G4ElectroMagneticField *m_inner; }; diff --git a/core/opengate_core/opengate_lib/GateField.h b/core/opengate_core/opengate_lib/GateField.h deleted file mode 100644 index 66b002c0e1..0000000000 --- a/core/opengate_core/opengate_lib/GateField.h +++ /dev/null @@ -1,61 +0,0 @@ -/* -------------------------------------------------- - Copyright (C): OpenGATE Collaboration - This software is distributed under the terms - of the GNU Lesser General Public Licence (LGPL) - See LICENSE.md for further details - -------------------------------------------------- */ - -#ifndef GateField_h -#define GateField_h - -#include "G4AffineTransform.hh" -#include "G4RotationMatrix.hh" -#include "G4ThreeVector.hh" -#include - -// Forward declaration to avoid including the full header. -class G4VSolid; - -// Helper base class for fields in GATE: -// - stores the world-to-local transforms of every physical -// placement of the logical volume the field is attached to. -// - provides the coordinate conversions needed to evaluate fields in each -// placement's local frame and rotate the result back to world. -class GateField { - -public: - // constructor — deltaChordMM is the chord-finder tolerance (delta_chord), - // used to compute the fallback-fatal distance internally. - GateField(const G4VSolid *solid, std::vector translations, - std::vector rotations, double deltaChordMM); - - // update the transforms (e.g. after a geometry change between runs) - void SetTransforms(std::vector translations, - std::vector rotations); - -protected: - // find the local coordinates of worldPoint in the containing placement, and - // return the world-to-local transform of that placement in outTransform. - G4ThreeVector - findContainingPlacement(const G4ThreeVector &worldPoint, - const G4AffineTransform *&outTransform) const; - - // rotate a field vector from local to world coordinates using the given - // transform. - static G4ThreeVector rotateToWorld(const G4ThreeVector &localField, - const G4AffineTransform &transform); - - const G4VSolid - *m_solid; // solid of the logical volume the field is attached to - - double - m_fallbackFatalDistanceMM; // computed from deltaChordMM at construction - - std::vector m_transforms; // one per physical placement - - // Upper bound on the radius of curvature [mm] used to derive - // the fallback-fatal distance. - static constexpr double kMaxCurvatureRadiusMM = 100.0e3; -}; - -#endif // GateField_h diff --git a/core/opengate_core/opengate_lib/GateField.cpp b/core/opengate_core/opengate_lib/GateFieldBase.cpp similarity index 55% rename from core/opengate_core/opengate_lib/GateField.cpp rename to core/opengate_core/opengate_lib/GateFieldBase.cpp index 7a1c4ec572..f4893fe0d4 100644 --- a/core/opengate_core/opengate_lib/GateField.cpp +++ b/core/opengate_core/opengate_lib/GateFieldBase.cpp @@ -5,36 +5,31 @@ See LICENSE.md for further details -------------------------------------------------- */ -#include "GateField.h" +#include "GateFieldBase.h" #include "G4VSolid.hh" -#include "globals.hh" // G4Exception +#include "globals.hh" #include -#include #include #include #include -#include // constructor -GateField::GateField(const G4VSolid *solid, - std::vector translations, - std::vector rotations, - double deltaChordMM) +GateFieldBase::GateFieldBase(const G4VSolid *solid, + std::vector translations, + std::vector rotations, + double deltaChordMM) : m_solid(solid), m_fallbackFatalDistanceMM( 5.0 * std::sqrt(8.0 * kMaxCurvatureRadiusMM * deltaChordMM)) { - // sanity-check the inputs before caching if (solid == nullptr) - throw std::invalid_argument("GateField: solid must not be null"); + throw std::invalid_argument("GateFieldBase: solid must not be null"); if (translations.size() != rotations.size() || translations.empty()) - throw std::invalid_argument("GateField: translations and rotations must be " - "non-empty and have the same size"); + throw std::invalid_argument("GateFieldBase: translations and rotations " + "must be non-empty and have the same size"); - // initial cache the world-to-local transforms for every physical placement of - // the physical volume m_transforms.reserve(translations.size()); for (std::size_t i = 0; i < translations.size(); ++i) m_transforms.emplace_back(rotations[i].inverse(), translations[i]); @@ -42,17 +37,13 @@ GateField::GateField(const G4VSolid *solid, // recache the world-to-local transforms (e.g. after a geometry change between // runs) -void GateField::SetTransforms(std::vector translations, - std::vector rotations) { - - // sanity-check the inputs before caching +void GateFieldBase::SetTransforms(std::vector translations, + std::vector rotations) { if (translations.size() != rotations.size() || translations.empty()) throw std::invalid_argument( - "GateField::SetTransforms: translations and rotations must be " + "GateFieldBase::SetTransforms: translations and rotations must be " "non-empty and have the same size"); - // recache the world-to-local transforms for every physical placement of the - // logical volume m_transforms.clear(); m_transforms.reserve(translations.size()); for (std::size_t i = 0; i < translations.size(); ++i) @@ -60,15 +51,10 @@ void GateField::SetTransforms(std::vector translations, } // find the local coordinates of worldPoint in the containing placement of the -// field's logical volume, and return the transform of that placement. -G4ThreeVector GateField::findContainingPlacement( +// field's logical volume +G4ThreeVector GateFieldBase::findContainingPlacement( const G4ThreeVector &worldPoint, const G4AffineTransform *&outTransform) const { - - // Loop over all placements once: - // - return immediately if the point is inside any placement; - // - otherwise accumulate the closest surface distance for the fallback. - // This avoids re-transforming the point in a second pass. std::size_t closestIdx = 0; double minDistToSurface = std::numeric_limits::infinity(); G4ThreeVector closestLocal{}; @@ -82,8 +68,6 @@ G4ThreeVector GateField::findContainingPlacement( return localPoint; } - // Fallback: track the placement whose surface is nearest. - // DistanceToIn is valid here because Inside() just returned kOutside. const double d = m_solid->DistanceToIn(localPoint); if (d < minDistToSurface) { minDistToSurface = d; @@ -92,34 +76,36 @@ G4ThreeVector GateField::findContainingPlacement( } } - // sanity check: if the closest surface is still too far, this is likely a - // real bug + // if the point is outside every placement, check if it's within the + // fallback-fatal distance which accounts for any reasonable field integrator + // overshoot due to tolerances. if not, this is very likely a bug in the + // geometry or field setup -> fatal if (minDistToSurface > m_fallbackFatalDistanceMM) { std::ostringstream msg; - msg << "GateField::findContainingPlacement: world point (" << worldPoint.x() - << ", " << worldPoint.y() << ", " << worldPoint.z() << ") mm is " - << minDistToSurface + msg << "GateFieldBase::findContainingPlacement: world point (" + << worldPoint.x() << ", " << worldPoint.y() << ", " << worldPoint.z() + << ") mm is " << minDistToSurface << " mm outside every cached placement of the field's solid — " << "well beyond any chord-finder overshoot.\n" << " Closest placement: index " << closestIdx << " (local point " << closestLocal.x() << ", " << closestLocal.y() << ", " << closestLocal.z() << ").\n" - << " Maximum allowed distance before fatal: " + << " Maximum allowed distance before fatal: " << m_fallbackFatalDistanceMM << " mm.\n" - << " This likely indicates a real bug in the geometry or field setup\n"; + << " This likely indicates a real bug in the geometry or field " + "setup.\n"; - G4Exception("GateField::findContainingPlacement", "GateField0001", + G4Exception("GateFieldBase::findContainingPlacement", "GateField0001", FatalException, msg.str().c_str()); } outTransform = &m_transforms[closestIdx]; - return closestLocal; } // rotate a field vector from local to world coordinates using the given // transform -G4ThreeVector GateField::rotateToWorld(const G4ThreeVector &localField, - const G4AffineTransform &transform) { +G4ThreeVector GateFieldBase::rotateToWorld(const G4ThreeVector &localField, + const G4AffineTransform &transform) { return transform.TransformAxis(localField); } diff --git a/core/opengate_core/opengate_lib/GateFieldBase.h b/core/opengate_core/opengate_lib/GateFieldBase.h new file mode 100644 index 0000000000..fa1c3adde3 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateFieldBase.h @@ -0,0 +1,81 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateFieldBase_h +#define GateFieldBase_h + +#include "G4AffineTransform.hh" +#include "G4RotationMatrix.hh" +#include "G4ThreeVector.hh" +#include "G4Types.hh" +#include + +class G4VSolid; + +// Shared base class for all GATE field types +// - stores the world-to-local transforms of every physical placement of the +// logical volume +// - handles the full coordinate transform journey for GetFieldValue +class GateFieldBase { +public: + // constructor + GateFieldBase(const G4VSolid *solid, std::vector translations, + std::vector rotations, double deltaChordMM); + + // update the transforms (e.g. after a geometry change between runs) + void SetTransforms(std::vector translations, + std::vector rotations); + +protected: + // LocalFieldFn is expected to be something like: + // void getLocalField(const G4double localPos[4], G4double *field); + // which fills the field array with the local field value at the given + // localPos. + template + + // Get the field value at the given world point by transforming to local + // coordinates, calling getLocalField to fill the field value in local frame, + // and then rotating the field vector(s) back to world frame + void applyTransforms(const G4double Point[4], G4double *field, + int nComponents, LocalFieldFn getLocalField) const { + const G4ThreeVector worldPt(Point[0], Point[1], Point[2]); + const G4AffineTransform *tr = nullptr; + const G4ThreeVector localPt = findContainingPlacement(worldPt, tr); + const G4double localPos[4] = {localPt.x(), localPt.y(), localPt.z(), + Point[3]}; + + getLocalField(localPos, field); + + for (int i = 0; i < nComponents; i += 3) { + const G4ThreeVector v = rotateToWorld( + G4ThreeVector(field[i], field[i + 1], field[i + 2]), *tr); + field[i] = v.x(); + field[i + 1] = v.y(); + field[i + 2] = v.z(); + } + } + +private: + // find the local coordinates of worldPoint in the containing placement of the + // field's logical volume + G4ThreeVector + findContainingPlacement(const G4ThreeVector &worldPoint, + const G4AffineTransform *&outTransform) const; + + // rotate a field vector from local to world coordinates using the given + // transform + static G4ThreeVector rotateToWorld(const G4ThreeVector &localField, + const G4AffineTransform &transform); + + const G4VSolid *m_solid; + double m_fallbackFatalDistanceMM; + std::vector m_transforms; + + static constexpr double kMaxCurvatureRadiusMM = 100.0e3; +}; + +#endif // GateFieldBase_h diff --git a/core/opengate_core/opengate_lib/GateMagneticField.cpp b/core/opengate_core/opengate_lib/GateMagneticField.cpp index d35f0035a0..d2e38e71e9 100644 --- a/core/opengate_core/opengate_lib/GateMagneticField.cpp +++ b/core/opengate_core/opengate_lib/GateMagneticField.cpp @@ -13,36 +13,17 @@ GateMagneticField::GateMagneticField(G4MagneticField *inner, std::vector translations, std::vector rotations, double deltaChordMM) - : G4MagneticField(), GateField(solid, std::move(translations), - std::move(rotations), deltaChordMM), + : G4MagneticField(), GateFieldBase(solid, std::move(translations), + std::move(rotations), deltaChordMM), m_inner(inner) {} -// override of G4MagneticField interface void GateMagneticField::GetFieldValue(const G4double Point[4], G4double *Bfield) const { + // localFieldFunc is a lambda that captures 'this' and calls + // m_inner->GetFieldValue + auto localFieldFunc = [this](const G4double pos[4], G4double *f) { + m_inner->GetFieldValue(pos, f); + }; - // get the global coordinates of the point and the transform of the containing - // placement - const G4ThreeVector worldPoint(Point[0], Point[1], Point[2]); - - // get the local coordinates of the point in the containing placement, and the - // transform of that placement - const G4AffineTransform *transform = nullptr; - const G4ThreeVector localPoint = - findContainingPlacement(worldPoint, transform); - const G4double localPos[4] = {localPoint.x(), localPoint.y(), localPoint.z(), - Point[3]}; - - // get the local field value from the inner field object - G4double localB[3] = {0.0, 0.0, 0.0}; - m_inner->GetFieldValue(localPos, localB); - - // rotate the field back to world coordinates - const G4ThreeVector worldB = - rotateToWorld({localB[0], localB[1], localB[2]}, *transform); - - // copy the result to the output array - Bfield[0] = worldB.x(); - Bfield[1] = worldB.y(); - Bfield[2] = worldB.z(); + applyTransforms(Point, Bfield, 3, localFieldFunc); } diff --git a/core/opengate_core/opengate_lib/GateMagneticField.h b/core/opengate_core/opengate_lib/GateMagneticField.h index 53bff33896..e05b3432b1 100644 --- a/core/opengate_core/opengate_lib/GateMagneticField.h +++ b/core/opengate_core/opengate_lib/GateMagneticField.h @@ -9,44 +9,26 @@ #define GateMagneticField_h #include "G4MagneticField.hh" -#include "G4RotationMatrix.hh" -#include "G4ThreeVector.hh" -#include "GateField.h" +#include "GateFieldBase.h" #include class G4VSolid; -// GATE wrapper for G4MagneticField that transforms the query point to the local -// coordinates of the physical volume(s) the field is attached to, then -// delegates to the inner field class to get the field value, and finally -// transforms the field vector back to world coordinates. - -// wrapper for an internal G4MagneticField converted to local coordinates -class GateMagneticField : public G4MagneticField, protected GateField { +// GATE wrapper for G4MagneticField +class GateMagneticField : public G4MagneticField, public GateFieldBase { public: - // constructor — deltaChordMM is the chord-finder tolerance (delta_chord), - // forwarded to GateField to compute the fallback-fatal distance internally. - GateMagneticField(G4MagneticField *inner, // wrapped inner field - const G4VSolid *solid, // solid of the logical volume the - // field is attached to - std::vector - translations, // translations of the physical placements - std::vector - rotations, // rotations of the physical placements of - // the logical volume + // constructor + GateMagneticField(G4MagneticField *inner, const G4VSolid *solid, + std::vector translations, + std::vector rotations, double deltaChordMM); - // override of G4MagneticField interface + // override GetFieldValue to apply the coordinate transforms void GetFieldValue(const G4double Point[4], G4double *Bfield) const override; - // forward override of GateField::SetTransforms - inline void SetTransforms(std::vector translations, - std::vector rotations) { - GateField::SetTransforms(std::move(translations), std::move(rotations)); - } - private: - G4MagneticField *m_inner; // wrapped inner field + // the inner field in the volume's local frame + G4MagneticField *m_inner; }; #endif // GateMagneticField_h From 22ad8ec620fcfb202c2f591347df7b8995976e41 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Wed, 6 May 2026 16:14:48 +0200 Subject: [PATCH 12/28] Implement grid-based mapped electric and magnetic fields Co-authored-by: Copilot --- core/opengate_core/opengate_core.cpp | 11 +- .../opengate_lib/GateGridInterpolator.cpp | 115 +++++++ .../opengate_lib/GateGridInterpolator.h | 59 ++++ .../opengate_lib/GateMappedElectricField.cpp | 26 ++ .../opengate_lib/GateMappedElectricField.h | 34 ++ .../opengate_lib/GateMappedFieldBase.cpp | 34 ++ .../opengate_lib/GateMappedFieldBase.h | 38 +++ .../opengate_lib/GateMappedMagneticField.cpp | 24 ++ .../opengate_lib/GateMappedMagneticField.h | 31 ++ .../opengate_lib/pyGateGridInterpolator.cpp | 21 ++ .../pyGateMappedElectricField.cpp | 59 ++++ .../pyGateMappedMagneticField.cpp | 59 ++++ opengate/geometry/fields.py | 313 ++++++++++-------- 13 files changed, 691 insertions(+), 133 deletions(-) create mode 100644 core/opengate_core/opengate_lib/GateGridInterpolator.cpp create mode 100644 core/opengate_core/opengate_lib/GateGridInterpolator.h create mode 100644 core/opengate_core/opengate_lib/GateMappedElectricField.cpp create mode 100644 core/opengate_core/opengate_lib/GateMappedElectricField.h create mode 100644 core/opengate_core/opengate_lib/GateMappedFieldBase.cpp create mode 100644 core/opengate_core/opengate_lib/GateMappedFieldBase.h create mode 100644 core/opengate_core/opengate_lib/GateMappedMagneticField.cpp create mode 100644 core/opengate_core/opengate_lib/GateMappedMagneticField.h create mode 100644 core/opengate_core/opengate_lib/pyGateGridInterpolator.cpp create mode 100644 core/opengate_core/opengate_lib/pyGateMappedElectricField.cpp create mode 100644 core/opengate_core/opengate_lib/pyGateMappedMagneticField.cpp diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index 29b274f948..d5ed056d10 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -215,8 +215,10 @@ void init_G4ChordFinder(py::module &); void init_GateMagneticField(py::module &); void init_GateElectroMagneticField(py::module &); void init_GateUniformElectroMagneticField(py::module &); - -// void init_GateMappedMagneticField(py::module &); +void init_GateGridInterpolator(py::module &); +void init_GateMappedMagneticField(py::module &); +void init_GateMappedElectricField(py::module &); +// void init_GateMappedElectroMagneticField(py::module &); // geometry/solids void init_G4Box(py::module &); @@ -596,7 +598,10 @@ PYBIND11_MODULE(opengate_core, m) { init_GateMagneticField(m); init_GateElectroMagneticField(m); init_GateUniformElectroMagneticField(m); - // init_GateMappedMagneticField(m); + init_GateGridInterpolator(m); + init_GateMappedMagneticField(m); + init_GateMappedElectricField(m); + // init_GateMappedElectroMagneticField(m); init_G4Box(m); init_G4Ellipsoid(m); diff --git a/core/opengate_core/opengate_lib/GateGridInterpolator.cpp b/core/opengate_core/opengate_lib/GateGridInterpolator.cpp new file mode 100644 index 0000000000..9f072fd09f --- /dev/null +++ b/core/opengate_core/opengate_lib/GateGridInterpolator.cpp @@ -0,0 +1,115 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include "GateGridInterpolator.h" + +#include +#include +#include + +// Constructor: validate that field values have the correct size +GateGridInterpolator::GateGridInterpolator(GridDefinition gridDef, + FieldValues fieldValues, + InterpolationMethod interpMethod) + : m_gridDef(gridDef), m_fieldValues(std::move(fieldValues)), + m_interpMethod(interpMethod) { + // Validate the grid definition and field values + const int N = gridDef.nx * gridDef.ny * gridDef.nz; + if ((int)m_fieldValues.Fx.size() != N || (int)m_fieldValues.Fy.size() != N || + (int)m_fieldValues.Fz.size() != N) { + throw std::invalid_argument("GateGridInterpolator: Field value arrays must " + "each have nx*ny*nz elements"); + } +} + +// interpolate the field value at a point given in local coordinates +void GateGridInterpolator::Interpolate(double x, double y, double z, + double *field) const { + const auto &g = m_gridDef; + double fx = (x - g.x0) / g.dx; + double fy = (y - g.y0) / g.dy; + double fz = (z - g.z0) / g.dz; + + switch (m_interpMethod) { + case InterpolationMethod::Nearest: + nearest(fx, fy, fz, field); + break; + case InterpolationMethod::Trilinear: + trilinear(fx, fy, fz, field); + break; + default: + throw std::runtime_error( + "GateGridInterpolator: Unknown interpolation method"); + } +} + +// Trilinear interpolation +void GateGridInterpolator::trilinear(double fx, double fy, double fz, + double *field) const { + // Compute the surrounding grid point indices + const int ix0 = std::clamp((int)std::floor(fx), 0, m_gridDef.nx - 1); + const int iy0 = std::clamp((int)std::floor(fy), 0, m_gridDef.ny - 1); + const int iz0 = std::clamp((int)std::floor(fz), 0, m_gridDef.nz - 1); + + const int ix1 = std::clamp(ix0 + 1, 0, m_gridDef.nx - 1); + const int iy1 = std::clamp(iy0 + 1, 0, m_gridDef.ny - 1); + const int iz1 = std::clamp(iz0 + 1, 0, m_gridDef.nz - 1); + + // Compute the fractional distances along each axis + const double tx = std::clamp(fx - ix0, 0.0, 1.0); + const double ty = std::clamp(fy - iy0, 0.0, 1.0); + const double tz = std::clamp(fz - iz0, 0.0, 1.0); + + double fieldX = + m_fieldValues.Fx[idx(ix0, iy0, iz0)] * (1 - tx) * (1 - ty) * (1 - tz) + + m_fieldValues.Fx[idx(ix1, iy0, iz0)] * tx * (1 - ty) * (1 - tz) + + m_fieldValues.Fx[idx(ix0, iy1, iz0)] * (1 - tx) * ty * (1 - tz) + + m_fieldValues.Fx[idx(ix1, iy1, iz0)] * tx * ty * (1 - tz) + + m_fieldValues.Fx[idx(ix0, iy0, iz1)] * (1 - tx) * (1 - ty) * tz + + m_fieldValues.Fx[idx(ix1, iy0, iz1)] * tx * (1 - ty) * tz + + m_fieldValues.Fx[idx(ix0, iy1, iz1)] * (1 - tx) * ty * tz + + m_fieldValues.Fx[idx(ix1, iy1, iz1)] * tx * ty * tz; + + double fieldY = + m_fieldValues.Fy[idx(ix0, iy0, iz0)] * (1 - tx) * (1 - ty) * (1 - tz) + + m_fieldValues.Fy[idx(ix1, iy0, iz0)] * tx * (1 - ty) * (1 - tz) + + m_fieldValues.Fy[idx(ix0, iy1, iz0)] * (1 - tx) * ty * (1 - tz) + + m_fieldValues.Fy[idx(ix1, iy1, iz0)] * tx * ty * (1 - tz) + + m_fieldValues.Fy[idx(ix0, iy0, iz1)] * (1 - tx) * (1 - ty) * tz + + m_fieldValues.Fy[idx(ix1, iy0, iz1)] * tx * (1 - ty) * tz + + m_fieldValues.Fy[idx(ix0, iy1, iz1)] * (1 - tx) * ty * tz + + m_fieldValues.Fy[idx(ix1, iy1, iz1)] * tx * ty * tz; + + double fieldZ = + m_fieldValues.Fz[idx(ix0, iy0, iz0)] * (1 - tx) * (1 - ty) * (1 - tz) + + m_fieldValues.Fz[idx(ix1, iy0, iz0)] * tx * (1 - ty) * (1 - tz) + + m_fieldValues.Fz[idx(ix0, iy1, iz0)] * (1 - tx) * ty * (1 - tz) + + m_fieldValues.Fz[idx(ix1, iy1, iz0)] * tx * ty * (1 - tz) + + m_fieldValues.Fz[idx(ix0, iy0, iz1)] * (1 - tx) * (1 - ty) * tz + + m_fieldValues.Fz[idx(ix1, iy0, iz1)] * tx * (1 - ty) * tz + + m_fieldValues.Fz[idx(ix0, iy1, iz1)] * (1 - tx) * ty * tz + + m_fieldValues.Fz[idx(ix1, iy1, iz1)] * tx * ty * tz; + + field[0] = fieldX; + field[1] = fieldY; + field[2] = fieldZ; +} + +// Nearest-neighbor interpolation +void GateGridInterpolator::nearest(double fx, double fy, double fz, + double *field) const { + // Compute the nearest grid point indices + const int ix = std::clamp((int)std::round(fx), 0, m_gridDef.nx - 1); + const int iy = std::clamp((int)std::round(fy), 0, m_gridDef.ny - 1); + const int iz = std::clamp((int)std::round(fz), 0, m_gridDef.nz - 1); + + const int i = idx(ix, iy, iz); + + field[0] = m_fieldValues.Fx[i]; + field[1] = m_fieldValues.Fy[i]; + field[2] = m_fieldValues.Fz[i]; +} diff --git a/core/opengate_core/opengate_lib/GateGridInterpolator.h b/core/opengate_core/opengate_lib/GateGridInterpolator.h new file mode 100644 index 0000000000..20e5791034 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateGridInterpolator.h @@ -0,0 +1,59 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateGridInterpolator_h +#define GateGridInterpolator_h + +#include + +// Utility class for interpolating field values defined on a regular 3D grid +class GateGridInterpolator { +public: + // container for grid definition parameters + class GridDefinition { + public: + int nx, ny, nz; // number of grid points along each axis + double x0, y0, z0; // origin of the grid in local coordinates + double dx, dy, dz; // grid spacing along each axis + }; + + // container for field values defined on the grid + // must be defined in lexicographical order x->y->z (x slowest, z fastest) + class FieldValues { + public: + std::vector Fx; // Field_x values + std::vector Fy; // Field_y values + std::vector Fz; // Field_z values + }; + + // interpolation methods + enum class InterpolationMethod { Nearest, Trilinear }; + + // constructor + GateGridInterpolator(GridDefinition gridDef, FieldValues fieldValues, + InterpolationMethod interpMethod); + + // interpolate the field value at a point given in local coordinates + void Interpolate(double x, double y, double z, double *field) const; + +private: + GridDefinition m_gridDef; + FieldValues m_fieldValues; + InterpolationMethod m_interpMethod; + + // compute the 1D index in the field arrays for given grid indices according + // to lexicographical order x->y->z (x slowest, z fastest) + inline int idx(int ix, int iy, int iz) const { + return ix * (m_gridDef.ny * m_gridDef.nz) + iy * m_gridDef.nz + iz; + } + + // interpolation methods + void trilinear(double fx, double fy, double fz, double *field) const; + void nearest(double fx, double fy, double fz, double *field) const; +}; + +#endif // GateGridInterpolator_h diff --git a/core/opengate_core/opengate_lib/GateMappedElectricField.cpp b/core/opengate_core/opengate_lib/GateMappedElectricField.cpp new file mode 100644 index 0000000000..05a0616a76 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateMappedElectricField.cpp @@ -0,0 +1,26 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include "GateMappedElectricField.h" + +GateMappedElectricField::GateMappedElectricField( + const G4VSolid *solid, std::vector translations, + std::vector rotations, double deltaChordMM, + GateGridInterpolator::GridDefinition gridDef, + GateGridInterpolator::FieldValues fieldValues, + GateGridInterpolator::InterpolationMethod interpMethod) + : G4ElectroMagneticField(), + GateMappedFieldBase(solid, std::move(translations), std::move(rotations), + deltaChordMM, gridDef, std::move(fieldValues), + interpMethod) {} + +void GateMappedElectricField::GetFieldValue(const G4double Point[4], + G4double *BEfield) const { + BEfield[0] = BEfield[1] = BEfield[2] = 0.0; + computeFieldValue(Point, BEfield + 3, + 3); // +3 points the output pointer to the E components +} diff --git a/core/opengate_core/opengate_lib/GateMappedElectricField.h b/core/opengate_core/opengate_lib/GateMappedElectricField.h new file mode 100644 index 0000000000..ea57ec775f --- /dev/null +++ b/core/opengate_core/opengate_lib/GateMappedElectricField.h @@ -0,0 +1,34 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateMappedElectricField_h +#define GateMappedElectricField_h + +#include "G4ElectroMagneticField.hh" +#include "GateMappedFieldBase.h" +#include + +class G4VSolid; + +// grid-based mapped electric field +class GateMappedElectricField : public G4ElectroMagneticField, + public GateMappedFieldBase { +public: + GateMappedElectricField( + const G4VSolid *solid, std::vector translations, + std::vector rotations, double deltaChordMM, + GateGridInterpolator::GridDefinition gridDef, + GateGridInterpolator::FieldValues fieldValues, + GateGridInterpolator::InterpolationMethod interpMethod); + + // returns [0, 0, 0, Ex, Ey, Ez] + void GetFieldValue(const G4double Point[4], G4double *BEfield) const override; + + G4bool DoesFieldChangeEnergy() const override { return true; } +}; + +#endif // GateMappedElectricField_h diff --git a/core/opengate_core/opengate_lib/GateMappedFieldBase.cpp b/core/opengate_core/opengate_lib/GateMappedFieldBase.cpp new file mode 100644 index 0000000000..457d3649c2 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateMappedFieldBase.cpp @@ -0,0 +1,34 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include "GateMappedFieldBase.h" + +// constructor +GateMappedFieldBase::GateMappedFieldBase( + const G4VSolid *solid, std::vector translations, + std::vector rotations, double deltaChordMM, + GateGridInterpolator::GridDefinition gridDef, + GateGridInterpolator::FieldValues fieldValues, + GateGridInterpolator::InterpolationMethod interpMethod) + : GateFieldBase(solid, std::move(translations), std::move(rotations), + deltaChordMM), + m_interpolator(gridDef, std::move(fieldValues), interpMethod) {} + +// Compute field value at world point: transform to local, interpolate, rotate +// back +void GateMappedFieldBase::computeFieldValue(const G4double Point[4], + G4double *field, + int nComponents) const { + + // localFieldFunc is a lambda that captures 'this' and calls + // m_interpolator.Interpolate with the local point coordinates + auto localFieldFunc = [this](const G4double pos[4], G4double *f) { + m_interpolator.Interpolate(pos[0], pos[1], pos[2], f); + }; + + applyTransforms(Point, field, nComponents, localFieldFunc); +} diff --git a/core/opengate_core/opengate_lib/GateMappedFieldBase.h b/core/opengate_core/opengate_lib/GateMappedFieldBase.h new file mode 100644 index 0000000000..8e286539e6 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateMappedFieldBase.h @@ -0,0 +1,38 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateMappedFieldBase_h +#define GateMappedFieldBase_h + +#include "GateFieldBase.h" +#include "GateGridInterpolator.h" +#include + +class G4VSolid; + +// base class for mapped fields that use grid interpolation +class GateMappedFieldBase : public GateFieldBase { +public: + // constructor + GateMappedFieldBase(const G4VSolid *solid, + std::vector translations, + std::vector rotations, + double deltaChordMM, + GateGridInterpolator::GridDefinition gridDef, + GateGridInterpolator::FieldValues fieldValues, + GateGridInterpolator::InterpolationMethod interpMethod); + +protected: + // transform Point to local frame, interpolate, rotate result back to world. + void computeFieldValue(const G4double Point[4], G4double *field, + int nComponents) const; + +private: + GateGridInterpolator m_interpolator; +}; + +#endif // GateMappedFieldBase_h diff --git a/core/opengate_core/opengate_lib/GateMappedMagneticField.cpp b/core/opengate_core/opengate_lib/GateMappedMagneticField.cpp new file mode 100644 index 0000000000..fc2d7a95d4 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateMappedMagneticField.cpp @@ -0,0 +1,24 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include "GateMappedMagneticField.h" + +GateMappedMagneticField::GateMappedMagneticField( + const G4VSolid *solid, std::vector translations, + std::vector rotations, double deltaChordMM, + GateGridInterpolator::GridDefinition gridDef, + GateGridInterpolator::FieldValues fieldValues, + GateGridInterpolator::InterpolationMethod interpMethod) + : G4MagneticField(), + GateMappedFieldBase(solid, std::move(translations), std::move(rotations), + deltaChordMM, gridDef, std::move(fieldValues), + interpMethod) {} + +void GateMappedMagneticField::GetFieldValue(const G4double Point[4], + G4double *Bfield) const { + computeFieldValue(Point, Bfield, 3); +} diff --git a/core/opengate_core/opengate_lib/GateMappedMagneticField.h b/core/opengate_core/opengate_lib/GateMappedMagneticField.h new file mode 100644 index 0000000000..dd537b3dd0 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateMappedMagneticField.h @@ -0,0 +1,31 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateMappedMagneticField_h +#define GateMappedMagneticField_h + +#include "G4MagneticField.hh" +#include "GateMappedFieldBase.h" +#include + +class G4VSolid; + +// grid-based mapped magnetic field. +class GateMappedMagneticField : public G4MagneticField, + public GateMappedFieldBase { +public: + GateMappedMagneticField( + const G4VSolid *solid, std::vector translations, + std::vector rotations, double deltaChordMM, + GateGridInterpolator::GridDefinition gridDef, + GateGridInterpolator::FieldValues fieldValues, + GateGridInterpolator::InterpolationMethod interpMethod); + + void GetFieldValue(const G4double Point[4], G4double *Bfield) const override; +}; + +#endif // GateMappedMagneticField_h diff --git a/core/opengate_core/opengate_lib/pyGateGridInterpolator.cpp b/core/opengate_core/opengate_lib/pyGateGridInterpolator.cpp new file mode 100644 index 0000000000..1b1915c2e5 --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateGridInterpolator.cpp @@ -0,0 +1,21 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +#include "GateGridInterpolator.h" + +namespace py = pybind11; + +// exposes the interpolation method enum +void init_GateGridInterpolator(py::module &m) { + py::enum_( + m, "GateGridInterpolationMethod") + .value("Nearest", GateGridInterpolator::InterpolationMethod::Nearest) + .value("Trilinear", GateGridInterpolator::InterpolationMethod::Trilinear) + .export_values(); +} diff --git a/core/opengate_core/opengate_lib/pyGateMappedElectricField.cpp b/core/opengate_core/opengate_lib/pyGateMappedElectricField.cpp new file mode 100644 index 0000000000..c00165ab2e --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateMappedElectricField.cpp @@ -0,0 +1,59 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include +#include +#include + +#include "G4ElectroMagneticField.hh" +#include "G4VSolid.hh" +#include "GateGridInterpolator.h" +#include "GateMappedElectricField.h" + +namespace py = pybind11; + +namespace { +using DoubleArray = + py::array_t; + +std::vector to_vec(const DoubleArray &arr) { + auto buf = arr.request(); + return std::vector(static_cast(buf.ptr), + static_cast(buf.ptr) + buf.size); +} +} // namespace + +void init_GateMappedElectricField(py::module &m) { + + py::class_>( + m, "GateMappedElectricField") + + .def(py::init([](const G4VSolid *solid, + std::vector translations, + std::vector rotations, + double delta_chord_mm, int nx, int ny, int nz, double x0, + double y0, double z0, double dx, double dy, double dz, + DoubleArray Ex, DoubleArray Ey, DoubleArray Ez, + GateGridInterpolator::InterpolationMethod interp) { + GateGridInterpolator::GridDefinition gridDef{nx, ny, nz, x0, y0, + z0, dx, dy, dz}; + GateGridInterpolator::FieldValues fieldValues{ + to_vec(Ex), to_vec(Ey), to_vec(Ez)}; + return new GateMappedElectricField(solid, translations, rotations, + delta_chord_mm, gridDef, + fieldValues, interp); + }), + py::arg("solid"), py::arg("translations"), py::arg("rotations"), + py::arg("delta_chord_mm"), py::arg("nx"), py::arg("ny"), + py::arg("nz"), py::arg("x0"), py::arg("y0"), py::arg("z0"), + py::arg("dx"), py::arg("dy"), py::arg("dz"), py::arg("Ex"), + py::arg("Ey"), py::arg("Ez"), py::arg("interpolation")) + + .def("SetTransforms", &GateMappedElectricField::SetTransforms, + py::arg("translations"), py::arg("rotations")); +} diff --git a/core/opengate_core/opengate_lib/pyGateMappedMagneticField.cpp b/core/opengate_core/opengate_lib/pyGateMappedMagneticField.cpp new file mode 100644 index 0000000000..49a1285ac7 --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateMappedMagneticField.cpp @@ -0,0 +1,59 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include +#include +#include + +#include "G4MagneticField.hh" +#include "G4VSolid.hh" +#include "GateGridInterpolator.h" +#include "GateMappedMagneticField.h" + +namespace py = pybind11; + +namespace { +using DoubleArray = + py::array_t; + +std::vector to_vec(const DoubleArray &arr) { + auto buf = arr.request(); + return std::vector(static_cast(buf.ptr), + static_cast(buf.ptr) + buf.size); +} +} // namespace + +void init_GateMappedMagneticField(py::module &m) { + + py::class_>( + m, "GateMappedMagneticField") + + .def(py::init([](const G4VSolid *solid, + std::vector translations, + std::vector rotations, + double delta_chord_mm, int nx, int ny, int nz, double x0, + double y0, double z0, double dx, double dy, double dz, + DoubleArray Bx, DoubleArray By, DoubleArray Bz, + GateGridInterpolator::InterpolationMethod interp) { + GateGridInterpolator::GridDefinition gridDef{nx, ny, nz, x0, y0, + z0, dx, dy, dz}; + GateGridInterpolator::FieldValues fieldValues{ + to_vec(Bx), to_vec(By), to_vec(Bz)}; + return new GateMappedMagneticField(solid, translations, rotations, + delta_chord_mm, gridDef, + fieldValues, interp); + }), + py::arg("solid"), py::arg("translations"), py::arg("rotations"), + py::arg("delta_chord_mm"), py::arg("nx"), py::arg("ny"), + py::arg("nz"), py::arg("x0"), py::arg("y0"), py::arg("z0"), + py::arg("dx"), py::arg("dy"), py::arg("dz"), py::arg("Bx"), + py::arg("By"), py::arg("Bz"), py::arg("interpolation")) + + .def("SetTransforms", &GateMappedMagneticField::SetTransforms, + py::arg("translations"), py::arg("rotations")); +} diff --git a/opengate/geometry/fields.py b/opengate/geometry/fields.py index ddebb62188..35c278153e 100644 --- a/opengate/geometry/fields.py +++ b/opengate/geometry/fields.py @@ -511,134 +511,183 @@ def to_dictionary(self): ) -# class MappedMagneticField(MagneticField): -# """Magnetic field defined by values on a regular 3D Cartesian grid. - -# Field lookup is performed entirely in C++ (no Python callbacks at tracking -# time) using either trilinear or nearest-neighbour interpolation. Points -# outside the grid extent return zero field. -# """ - -# # ! TODO's (@srmarcballestero): -# # ! - refactor the _create_field() function to make it simpler. -# # ! - warning when using nearest interpolation with coarse grids -# # ! - implement a MappedElectroMagneticField, and MappedElectricField. -# # ! + for the electromagnetic case, separate grids should be used. -# # ! + a mother class should abstract the core functionality. - -# # hints for IDE -# field_matrix: np.ndarray -# interpolation: str - -# user_info_defaults = { -# "field_matrix": ( -# None, -# { -# "doc": ( -# "2D numpy array of shape (N, 6) with columns [x, y, z, Bx, By, Bz] " -# "on a regular Cartesian grid. Geant4 units should be used." -# ), -# }, -# ), -# "interpolation": ( -# "trilinear", -# { -# "doc": "Interpolation method: 'trilinear' (default) or 'nearest'.", -# }, -# ), -# } - -# def __init__(self, *args, **kwargs) -> None: -# super().__init__(*args, **kwargs) - -# def _create_field(self) -> None: -# if self.field_matrix is None: -# raise ValueError("field_matrix must be provided for MappedMagneticField") - -# mat = np.asarray(self.field_matrix, dtype=np.float64) -# if mat.ndim != 2 or mat.shape[1] != 6: -# raise ValueError( -# "field_matrix must be a 2D array with shape (N, 6) " -# "containing columns [x, y, z, Bx, By, Bz]" -# ) - -# # Separate coordinates and field values -# positions = mat[:, :3] -# B_values = mat[:, 3:] - -# # Sort in lexicographical order: x slowest, z fastest — matches C++ flat index. -# sort_idx = np.lexsort((positions[:, 2], positions[:, 1], positions[:, 0])) -# positions = positions[sort_idx] -# B_values = B_values[sort_idx] - -# # Round to suppress floating-point noise before uniqueness check. -# x_vals = np.unique(np.round(positions[:, 0], 10)) -# y_vals = np.unique(np.round(positions[:, 1], 10)) -# z_vals = np.unique(np.round(positions[:, 2], 10)) -# nx, ny, nz = len(x_vals), len(y_vals), len(z_vals) - -# # --- Grid validation --- -# if len(positions) != nx * ny * nz: -# raise ValueError( -# f"field_matrix does not define a complete regular 3D grid: " -# f"expected {nx}*{ny}*{nz}={nx * ny * nz} points, got {len(positions)}" -# ) -# if nx > 2 and not np.allclose(np.diff(x_vals), x_vals[1] - x_vals[0]): -# raise ValueError("field_matrix x-axis does not have uniform spacing") -# if ny > 2 and not np.allclose(np.diff(y_vals), y_vals[1] - y_vals[0]): -# raise ValueError("field_matrix y-axis does not have uniform spacing") -# if nz > 2 and not np.allclose(np.diff(z_vals), z_vals[1] - z_vals[0]): -# raise ValueError("field_matrix z-axis does not have uniform spacing") -# if nx < 2 or ny < 2 or nz < 2: -# raise ValueError( -# "field_matrix must have at least 2 unique points along each axis " -# "(no degenerate axes allowed)" -# ) -# # ---------------------- - -# x0, y0, z0 = float(x_vals[0]), float(y_vals[0]), float(z_vals[0]) -# dx = float(x_vals[1] - x_vals[0]) -# dy = float(y_vals[1] - y_vals[0]) -# dz = float(z_vals[1] - z_vals[0]) - -# interp_map = { -# "trilinear": g4.GateMappedMagneticFieldInterpolation.Trilinear, -# "nearest": g4.GateMappedMagneticFieldInterpolation.Nearest, -# } -# if self.interpolation not in interp_map: -# raise ValueError( -# f"Unknown interpolation '{self.interpolation}'. " -# f"Available options are 'trilinear' or 'nearest'." -# ) - -# # Collect one local-to-world transform per physical placement. -# translations_np, rotations_np = get_transform_world_to_local( -# self._field_volume_obj -# ) - -# g4_translations = [vec_np_as_g4(t) for t in translations_np] -# g4_rotations = [rot_np_as_g4(r) for r in rotations_np] - -# self.g4_field = g4.GateMappedMagneticField( -# B_values[:, 0], B_values[:, 1], B_values[:, 2], -# nx, ny, nz, -# x0, y0, z0, -# dx, dy, dz, -# interp_map[self.interpolation], -# g4_translations, g4_rotations, -# ) - -# # Serialization -# # def to_dictionary(self): -# # d = super().to_dictionary() -# # if self.field_matrix is not None: -# # d["field_matrix"] = np.asarray(self.field_matrix).tolist() -# # return d - -# # def from_dictionary(self, d): -# # super().from_dictionary(d) -# # if "field_matrix" in d and d["field_matrix"] is not None: -# # self.field_matrix = np.asarray(d["field_matrix"]) +def _parse_field_matrix(mat, class_name): + """Parse a field matrix into grid parameters and sorted field value arrays. + + mat must be a 2D array with columns [x, y, z, Fx, Fy, Fz] on a + regular Cartesian grid in Geant4 units, sorted in any order. + + Returns (nx, ny, nz, x0, y0, z0, dx, dy, dz, field_cols) where field_cols + is a list of n_field_cols 1D arrays in lexicographical x->y->z order. + """ + mat = np.asarray(mat, dtype=np.float64) + expected_cols = 6 + if mat.ndim != 2 or mat.shape[1] != expected_cols: + raise ValueError( + f"{class_name}: field_matrix must be a 2D array with {expected_cols} " + f"columns [x, y, z, Fx, Fy, Fz], " + f"got shape {mat.shape}" + ) + + positions = mat[:, :3] + field_values = mat[:, 3:] + + sort_idx = np.lexsort((positions[:, 2], positions[:, 1], positions[:, 0])) + positions = positions[sort_idx] + field_values = field_values[sort_idx] + + x_vals = np.unique(np.round(positions[:, 0], 10)) + y_vals = np.unique(np.round(positions[:, 1], 10)) + z_vals = np.unique(np.round(positions[:, 2], 10)) + nx, ny, nz = len(x_vals), len(y_vals), len(z_vals) + + if len(positions) != nx * ny * nz: + raise ValueError( + f"{class_name}: field_matrix does not define a complete regular 3D grid: " + f"expected {nx}*{ny}*{nz}={nx * ny * nz} points, got {len(positions)}" + ) + for axis, vals, label in ((x_vals, nx, "x"), (y_vals, ny, "y"), (z_vals, nz, "z")): + if vals < 2: + raise ValueError( + f"{class_name}: field_matrix must have at least 2 unique points " + f"along the {label}-axis" + ) + if vals > 2 and not np.allclose(np.diff(axis), axis[1] - axis[0]): + raise ValueError( + f"{class_name}: field_matrix {label}-axis does not have uniform spacing" + ) + + x0, y0, z0 = float(x_vals[0]), float(y_vals[0]), float(z_vals[0]) + dx = float(x_vals[1] - x_vals[0]) + dy = float(y_vals[1] - y_vals[0]) + dz = float(z_vals[1] - z_vals[0]) + + field_cols = [field_values[:, i] for i in range(3)] + return nx, ny, nz, x0, y0, z0, dx, dy, dz, field_cols + + +_interp_map = { + "trilinear": g4.GateGridInterpolationMethod.Trilinear, + "nearest": g4.GateGridInterpolationMethod.Nearest, +} + +_mapped_field_user_info = { + "field_matrix": ( + None, + { + "doc": ( + "2D numpy array on a regular Cartesian grid in Geant4 units. " + "Structure: [[x, y, z, field components...], ...]. " + ), + }, + ), + "interpolation": ( + "trilinear", + { + "doc": "Interpolation method: 'trilinear' (default) or 'nearest'.", + }, + ), +} + + +class MappedMagneticField(MagneticField): + """Magnetic field defined by values on a regular 3D Cartesian grid.""" + + field_matrix: np.ndarray + interpolation: str + + user_info_defaults = _mapped_field_user_info + + def create_field_manager(self, volume_obj) -> g4.G4FieldManager: + if self.field_matrix is None: + raise ValueError("field_matrix must be provided for MappedMagneticField") + if self.interpolation not in _interp_map: + raise ValueError( + f"Unknown interpolation '{self.interpolation}'. " + f"Choose 'trilinear' or 'nearest'." + ) + self._field_volume_obj = volume_obj + nx, ny, nz, x0, y0, z0, dx, dy, dz, (Bx, By, Bz) = _parse_field_matrix( + self.field_matrix, "MappedMagneticField" + ) + g4_translations, g4_rotations = self._make_g4_transforms() + gate_field = g4.GateMappedMagneticField( + self._field_volume_obj.g4_solid, + g4_translations, + g4_rotations, + self.delta_chord, + nx, + ny, + nz, + x0, + y0, + z0, + dx, + dy, + dz, + Bx, + By, + Bz, + _interp_map[self.interpolation], + ) + return self._build_field_manager( + None, gate_field, g4.G4Mag_UsualEqRhs, 6, volume_obj + ) + + +class MappedElectricField(ElectricField): + """Electric field defined by values on a regular 3D Cartesian grid.""" + + field_matrix: np.ndarray + interpolation: str + + user_info_defaults = _mapped_field_user_info + + def create_field_manager(self, volume_obj) -> g4.G4FieldManager: + if self.field_matrix is None: + raise ValueError("field_matrix must be provided for MappedElectricField") + if self.interpolation not in _interp_map: + raise ValueError( + f"Unknown interpolation '{self.interpolation}'. " + f"Choose 'trilinear' or 'nearest'." + ) + self._field_volume_obj = volume_obj + nx, ny, nz, x0, y0, z0, dx, dy, dz, (Ex, Ey, Ez) = _parse_field_matrix( + self.field_matrix, "MappedElectricField" + ) + g4_translations, g4_rotations = self._make_g4_transforms() + gate_field = g4.GateMappedElectricField( + self._field_volume_obj.g4_solid, + g4_translations, + g4_rotations, + self.delta_chord, + nx, + ny, + nz, + x0, + y0, + z0, + dx, + dy, + dz, + Ex, + Ey, + Ez, + _interp_map[self.interpolation], + ) + return self._build_field_manager( + None, gate_field, g4.G4EqMagElectricField, 8, volume_obj + ) + + +class MappedElectroMagneticField(ElectroMagneticField): + """Electromagnetic field defined by values on a regular 3D Cartesian grid. + + FIXME: implement this class. + """ + + def create_field_manager(self, volume_obj) -> g4.G4FieldManager: + raise NotImplementedError("MappedElectroMagneticField is not yet implemented.") field_types = { @@ -646,11 +695,13 @@ def to_dictionary(self): "QuadrupoleMagneticField": QuadrupoleMagneticField, "SextupoleMagneticField": SextupoleMagneticField, "CustomMagneticField": CustomMagneticField, - # "MappedMagneticField": MappedMagneticField, + "MappedMagneticField": MappedMagneticField, "UniformElectricField": UniformElectricField, "CustomElectricField": CustomElectricField, + "MappedElectricField": MappedElectricField, "UniformElectroMagneticField": UniformElectroMagneticField, "CustomElectroMagneticField": CustomElectroMagneticField, + "MappedElectroMagneticField": MappedElectroMagneticField, } process_cls(FieldBase) @@ -659,10 +710,12 @@ def to_dictionary(self): process_cls(QuadrupoleMagneticField) process_cls(SextupoleMagneticField) process_cls(CustomMagneticField) -# process_cls(MappedMagneticField) +process_cls(MappedMagneticField) process_cls(ElectroMagneticField) process_cls(ElectricField) process_cls(UniformElectricField) process_cls(CustomElectricField) +process_cls(MappedElectricField) process_cls(UniformElectroMagneticField) process_cls(CustomElectroMagneticField) +process_cls(MappedElectroMagneticField) From 60f893dafd73a361e7e3f7ef5ccac267c4d19d30 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Wed, 6 May 2026 16:14:48 +0200 Subject: [PATCH 13/28] Implement grid-based mapped electric and magnetic fields Co-authored-by: Copilot --- .../test099_fields_mapped_vs_uniform_B.py | 158 +++++++++++++++++ .../test099_fields_mapped_vs_uniform_E.py | 161 ++++++++++++++++++ 2 files changed, 319 insertions(+) create mode 100644 opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_B.py create mode 100644 opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py diff --git a/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_B.py b/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_B.py new file mode 100644 index 0000000000..a98668051f --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_B.py @@ -0,0 +1,158 @@ +#!/usr/bin/env python3 +""" +Test 099 — MappedMagneticField with uniform grid vs UniformMagneticField. + +Places two boxes side-by-side: one with UniformMagneticField, the other with +a MappedMagneticField whose grid contains a constant B value everywhere. A +proton source fires into each box. + +Checks that: + - Exit positions and energies agree within numerical tolerance. + - Both results agree with the analytical cyclotron radius. +""" + +import itertools + +import numpy as np +import uproot + +import opengate as gate +from opengate.geometry import fields +from opengate.tests import utility + +from test099_fields_helpers import ( + g4_m, + g4_cm, + g4_mm, + g4_tesla, + g4_MeV, + g4_eplus, + PROTON_MASS, + cyclotron_radius, +) + +VISU = False + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, output_folder="test099_fields") + + By = 5 * g4_tesla + T = 10 * g4_MeV + box_half = 250 * g4_mm # half-size of the 50 cm box + + # Build a uniform grid covering the box exactly. + # A 2x2x2 grid is the minimum for trilinear interpolation; on a uniform + # field it gives machine-precision exact results at every interior point. + corners = list(itertools.product([-box_half, box_half], repeat=3)) + field_matrix = np.array([[x, y, z, 0.0, By, 0.0] for x, y, z in corners]) + + sim = gate.Simulation() + sim.g4_verbose = False + sim.visu = False + sim.number_of_threads = 1 + sim.random_seed = 42 + sim.output_dir = paths.output + + if VISU: + sim.visu = True + sim.visu_type = "qt" + sim.visu_commands.append("/vis/scene/endOfEventAction accumulate") + sim.visu_commands.append("/vis/scene/add/trajectories smooth") + sim.visu_commands.append("/vis/scene/add/magneticField 15 fullArrow") + + world = sim.world + world.size = [3 * g4_m, 1 * g4_m, 1 * g4_m] + world.material = "G4_Galactic" + + # --- Box with UniformMagneticField (reference) --- + box_ref = sim.add_volume("BoxVolume", "box_ref") + box_ref.size = [50 * g4_cm, 50 * g4_cm, 50 * g4_cm] + box_ref.material = "G4_Galactic" + box_ref.translation = [-80 * g4_cm, 0, 0] + + uniform_field = fields.UniformMagneticField(name="B_uniform") + uniform_field.field_vector = [0, By, 0] + box_ref.add_field(uniform_field) + + # --- Box with MappedMagneticField (under test) --- + box_mapped = sim.add_volume("BoxVolume", "box_mapped") + box_mapped.size = [50 * g4_cm, 50 * g4_cm, 50 * g4_cm] + box_mapped.material = "G4_Galactic" + box_mapped.translation = [80 * g4_cm, 0, 0] + + mapped_field = fields.MappedMagneticField(name="B_mapped") + mapped_field.field_matrix = field_matrix + mapped_field.interpolation = "nearest" + box_mapped.add_field(mapped_field) + + # --- Sources --- + for name, translation in [ + ("src_ref", [-80 * g4_cm, 0, -100 * g4_cm]), + ("src_mapped", [80 * g4_cm, 0, -100 * g4_cm]), + ]: + src = sim.add_source("GenericSource", name) + src.particle = "proton" + src.n = 1 + src.energy.type = "mono" + src.energy.mono = T + src.position.type = "point" + src.position.translation = translation + src.direction.type = "momentum" + src.direction.momentum = [0, 0, 1] + + # --- Phase space actors --- + for name, box_name, filename in [ + ("phsp_ref", "box_ref", "test099_mapped_vs_uniform_ref.root"), + ("phsp_mapped", "box_mapped", "test099_mapped_vs_uniform_mapped.root"), + ]: + phsp = sim.add_actor("PhaseSpaceActor", name) + phsp.attached_to = box_name + phsp.attributes = ["PostKineticEnergy", "PostPosition"] + phsp.output_filename = paths.output / filename + phsp.steps_to_store = "exiting" + + sim.run() + + # --- Read results --- + df_ref = uproot.open(str(paths.output / "test099_mapped_vs_uniform_ref.root"))[ + "phsp_ref;1" + ].arrays(library="pd") + df_mapped = uproot.open( + str(paths.output / "test099_mapped_vs_uniform_mapped.root") + )["phsp_mapped;1"].arrays(library="pd") + + # Shift positions to each box's local frame for comparison + ref_x = df_ref["PostPosition_X"].values - (-80 * g4_cm) + mapped_x = df_mapped["PostPosition_X"].values - (80 * g4_cm) + ref_z = df_ref["PostPosition_Z"].values + mapped_z = df_mapped["PostPosition_Z"].values + ref_KE = df_ref["PostKineticEnergy"].values + mapped_KE = df_mapped["PostKineticEnergy"].values + + r_TOL = 1e-3 * g4_mm + e_TOL = 1e-3 * g4_MeV + + # Uniform vs mapped comparison + dx = np.abs(ref_x - mapped_x) + dz = np.abs(ref_z - mapped_z) + dKE = np.abs(ref_KE - mapped_KE) + + is_ok_x = np.all(dx < r_TOL) + is_ok_z = np.all(dz < r_TOL) + is_ok_e = np.all(dKE < e_TOL) + + print(f"Uniform vs mapped — dx max : {dx.max() / g4_mm:.6f} mm — OK: {is_ok_x}") + print(f"Uniform vs mapped — dz max : {dz.max() / g4_mm:.6f} mm — OK: {is_ok_z}") + print(f"Uniform vs mapped — dKE max: {dKE.max() / g4_MeV:.6f} MeV — OK: {is_ok_e}") + + # Analytical cyclotron check on the mapped result + r = cyclotron_radius(T, By, PROTON_MASS, 1 * g4_eplus) + cx, cz = 0.0 - r, -box_half + residual = np.sqrt((mapped_x - cx) ** 2 + (mapped_z - cz) ** 2) - r + is_ok_analytical = np.all(np.abs(residual) < r_TOL) + print( + f"Mapped vs analytical circle residual: {residual / g4_mm} mm — OK: {is_ok_analytical}" + ) + + is_ok = is_ok_x and is_ok_z and is_ok_e and is_ok_analytical + utility.test_ok(is_ok) diff --git a/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py b/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py new file mode 100644 index 0000000000..a0129bf815 --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py @@ -0,0 +1,161 @@ +#!/usr/bin/env python3 +""" +Test 099 — MappedMagneticField with uniform grid vs UniformMagneticField. + +Places two boxes side-by-side: one with UniformMagneticField, the other with +a MappedMagneticField whose grid contains a constant B value everywhere. A +proton source fires into each box. + +Checks that: + - Exit positions and energies agree within numerical tolerance. + - Both results agree with the analytical cyclotron radius. +""" + +import itertools + +import numpy as np +import uproot + +import opengate as gate +from opengate.geometry import fields +from opengate.tests import utility + +from test099_fields_helpers import ( + g4_m, + g4_cm, + g4_mm, + g4_volt, + g4_MeV, + g4_eplus, + PROTON_MASS, + cyclotron_radius, +) + +VISU = False + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, output_folder="test099_fields") + + Ex = 1e8 * g4_volt / g4_m + T = 10 * g4_MeV + q = 1 * g4_eplus + box_half = 250 * g4_mm # half-size of the 50 cm box + + # Build a uniform grid covering the box exactly. + # A 2x2x2 grid is the minimum for trilinear interpolation; on a uniform + # field it gives machine-precision exact results at every interior point. + corners = list(itertools.product([-1 * box_half, 1 * box_half], repeat=3)) + field_matrix = np.array([[x, y, z, Ex, 0.0, 0.0] for x, y, z in corners]) + + sim = gate.Simulation() + sim.g4_verbose = False + sim.visu = False + sim.number_of_threads = 1 + sim.random_seed = 42 + sim.output_dir = paths.output + + if VISU: + sim.visu = True + sim.visu_type = "qt" + sim.visu_commands.append("/vis/scene/endOfEventAction accumulate") + sim.visu_commands.append("/vis/scene/add/trajectories smooth") + sim.visu_commands.append("/vis/scene/add/electricField 15 fullArrow") + + world = sim.world + world.size = [3 * g4_m, 1 * g4_m, 1 * g4_m] + world.material = "G4_Galactic" + + # --- Box with UniformElectricField (reference) --- + box_ref = sim.add_volume("BoxVolume", "box_ref") + box_ref.size = [50 * g4_cm, 50 * g4_cm, 50 * g4_cm] + box_ref.material = "G4_Galactic" + box_ref.translation = [-80 * g4_cm, 0, 0] + + uniform_field = fields.UniformElectricField(name="E_uniform") + uniform_field.field_vector = [Ex, 0, 0] + box_ref.add_field(uniform_field) + + # --- Box with MappedElectricField (under test) --- + box_mapped = sim.add_volume("BoxVolume", "box_mapped") + box_mapped.size = [50 * g4_cm, 50 * g4_cm, 50 * g4_cm] + box_mapped.material = "G4_Galactic" + box_mapped.translation = [80 * g4_cm, 0, 0] + + mapped_field = fields.MappedElectricField(name="E_mapped") + mapped_field.field_matrix = field_matrix + mapped_field.interpolation = "nearest" + box_mapped.add_field(mapped_field) + + # --- Sources --- + for name, translation in [ + ("src_ref", [-80 * g4_cm, 0, -100 * g4_cm]), + ("src_mapped", [80 * g4_cm, 0, -100 * g4_cm]), + ]: + src = sim.add_source("GenericSource", name) + src.particle = "proton" + src.n = 1 + src.energy.type = "mono" + src.energy.mono = T + src.position.type = "point" + src.position.translation = translation + src.direction.type = "momentum" + src.direction.momentum = [0, 0, 1] + + # --- Phase space actors --- + for name, box_name, filename in [ + ("phsp_ref", "box_ref", "test099_mapped_vs_uniform_ref.root"), + ("phsp_mapped", "box_mapped", "test099_mapped_vs_uniform_mapped.root"), + ]: + phsp = sim.add_actor("PhaseSpaceActor", name) + phsp.attached_to = box_name + phsp.attributes = ["PostKineticEnergy", "PostPosition"] + phsp.output_filename = paths.output / filename + phsp.steps_to_store = "exiting" + + sim.run() + + # --- Read results --- + df_ref = uproot.open(str(paths.output / "test099_mapped_vs_uniform_ref.root"))[ + "phsp_ref;1" + ].arrays(library="pd") + df_mapped = uproot.open( + str(paths.output / "test099_mapped_vs_uniform_mapped.root") + )["phsp_mapped;1"].arrays(library="pd") + + # Shift positions to each box's local frame for comparison + ref_x = df_ref["PostPosition_X"].values - (-80 * g4_cm) + mapped_x = df_mapped["PostPosition_X"].values - (80 * g4_cm) + ref_y = df_ref["PostPosition_Y"].values + mapped_y = df_mapped["PostPosition_Y"].values + ref_z = df_ref["PostPosition_Z"].values + mapped_z = df_mapped["PostPosition_Z"].values + ref_KE = df_ref["PostKineticEnergy"].values + mapped_KE = df_mapped["PostKineticEnergy"].values + + r_TOL = 1e-3 * g4_mm + e_TOL = 1e-3 * g4_MeV + + # Uniform vs mapped comparison + dx = np.abs(ref_x - mapped_x) + dy = np.abs(ref_y - mapped_y) + dz = np.abs(ref_z - mapped_z) + dKE = np.abs(ref_KE - mapped_KE) + + is_ok_x = np.all(dx < r_TOL) + is_ok_y = np.all(dy < r_TOL) + is_ok_z = np.all(dz < r_TOL) + is_ok_e = np.all(dKE < e_TOL) + + print(f"Uniform vs mapped - dx max : {dx.max() / g4_mm:.6f} mm - OK: {is_ok_x}") + print(f"Uniform vs mapped - dy max : {dy.max() / g4_mm:.6f} mm - OK: {is_ok_y}") + print(f"Uniform vs mapped - dz max : {dz.max() / g4_mm:.6f} mm - OK: {is_ok_z}") + print(f"Uniform vs mapped - dKE max: {dKE.max() / g4_MeV:.6f} MeV - OK: {is_ok_e}") + + # Analytical check + x_from_energy = (mapped_KE - T) / (q * Ex) # ΔKE = qEΔx + is_ok_analytical = np.all(np.abs(mapped_x - x_from_energy) < r_TOL) + print( + f"Mapped vs analytical - dx max: {np.abs(mapped_x - x_from_energy).max() / g4_mm:.6f} mm - OK: {is_ok_analytical}" + ) + is_ok = is_ok_x and is_ok_z and is_ok_e and is_ok_analytical + utility.test_ok(is_ok) From aca9dad8e147bcc2d5829c356afdb78b3f090d69 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Wed, 6 May 2026 16:30:17 +0200 Subject: [PATCH 14/28] Implement mapped combined electric and magnetic fields Co-authored-by: Copilot --- core/opengate_core/opengate_core.cpp | 4 +- .../GateMappedElectroMagneticField.cpp | 41 +++++++ .../GateMappedElectroMagneticField.h | 41 +++++++ .../pyGateMappedElectroMagneticField.cpp | 70 +++++++++++ opengate/geometry/fields.py | 110 ++++++++++++++++-- 5 files changed, 255 insertions(+), 11 deletions(-) create mode 100644 core/opengate_core/opengate_lib/GateMappedElectroMagneticField.cpp create mode 100644 core/opengate_core/opengate_lib/GateMappedElectroMagneticField.h create mode 100644 core/opengate_core/opengate_lib/pyGateMappedElectroMagneticField.cpp diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index d5ed056d10..e9a6c1ec1a 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -218,7 +218,7 @@ void init_GateUniformElectroMagneticField(py::module &); void init_GateGridInterpolator(py::module &); void init_GateMappedMagneticField(py::module &); void init_GateMappedElectricField(py::module &); -// void init_GateMappedElectroMagneticField(py::module &); +void init_GateMappedElectroMagneticField(py::module &); // geometry/solids void init_G4Box(py::module &); @@ -601,7 +601,7 @@ PYBIND11_MODULE(opengate_core, m) { init_GateGridInterpolator(m); init_GateMappedMagneticField(m); init_GateMappedElectricField(m); - // init_GateMappedElectroMagneticField(m); + init_GateMappedElectroMagneticField(m); init_G4Box(m); init_G4Ellipsoid(m); diff --git a/core/opengate_core/opengate_lib/GateMappedElectroMagneticField.cpp b/core/opengate_core/opengate_lib/GateMappedElectroMagneticField.cpp new file mode 100644 index 0000000000..ad3e457135 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateMappedElectroMagneticField.cpp @@ -0,0 +1,41 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include "GateMappedElectroMagneticField.h" + +GateMappedElectroMagneticField::GateMappedElectroMagneticField( + const G4VSolid *solid, std::vector translations, + std::vector rotations, double deltaChordMM, + GateGridInterpolator::GridDefinition gridDefB, + GateGridInterpolator::FieldValues fieldValuesB, + GateGridInterpolator::GridDefinition gridDefE, + GateGridInterpolator::FieldValues fieldValuesE, + GateGridInterpolator::InterpolationMethod interpMethod) + : G4ElectroMagneticField(), + GateFieldBase(solid, std::move(translations), std::move(rotations), + deltaChordMM), + m_interpolator_B(gridDefB, std::move(fieldValuesB), interpMethod), + m_interpolator_E(gridDefE, std::move(fieldValuesE), interpMethod) {} + +void GateMappedElectroMagneticField::GetFieldValue(const G4double Point[4], + G4double *BEfield) const { + // localFieldFunc interpolates B and E separately, then combines them + auto localFieldFunc = [this](const G4double pos[4], G4double *f) { + double B[3] = {0.0, 0.0, 0.0}; + double E[3] = {0.0, 0.0, 0.0}; + m_interpolator_B.Interpolate(pos[0], pos[1], pos[2], B); + m_interpolator_E.Interpolate(pos[0], pos[1], pos[2], E); + f[0] = B[0]; + f[1] = B[1]; + f[2] = B[2]; + f[3] = E[0]; + f[4] = E[1]; + f[5] = E[2]; + }; + + applyTransforms(Point, BEfield, 6, localFieldFunc); +} diff --git a/core/opengate_core/opengate_lib/GateMappedElectroMagneticField.h b/core/opengate_core/opengate_lib/GateMappedElectroMagneticField.h new file mode 100644 index 0000000000..1a4cc0dbe8 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateMappedElectroMagneticField.h @@ -0,0 +1,41 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateMappedElectroMagneticField_h +#define GateMappedElectroMagneticField_h + +#include "G4ElectroMagneticField.hh" +#include "GateFieldBase.h" +#include "GateGridInterpolator.h" +#include + +class G4VSolid; + +// grid-based mapped electromagnetic field with separate B and E grids +class GateMappedElectroMagneticField : public G4ElectroMagneticField, + public GateFieldBase { +public: + GateMappedElectroMagneticField( + const G4VSolid *solid, std::vector translations, + std::vector rotations, double deltaChordMM, + GateGridInterpolator::GridDefinition gridDefB, + GateGridInterpolator::FieldValues fieldValuesB, + GateGridInterpolator::GridDefinition gridDefE, + GateGridInterpolator::FieldValues fieldValuesE, + GateGridInterpolator::InterpolationMethod interpMethod); + + // returns [Bx, By, Bz, Ex, Ey, Ez] + void GetFieldValue(const G4double Point[4], G4double *BEfield) const override; + + G4bool DoesFieldChangeEnergy() const override { return true; } + +private: + GateGridInterpolator m_interpolator_B; + GateGridInterpolator m_interpolator_E; +}; + +#endif // GateMappedElectroMagneticField_h \ No newline at end of file diff --git a/core/opengate_core/opengate_lib/pyGateMappedElectroMagneticField.cpp b/core/opengate_core/opengate_lib/pyGateMappedElectroMagneticField.cpp new file mode 100644 index 0000000000..6f2bb46d97 --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateMappedElectroMagneticField.cpp @@ -0,0 +1,70 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include +#include +#include + +#include "G4ElectroMagneticField.hh" +#include "G4VSolid.hh" +#include "GateGridInterpolator.h" +#include "GateMappedElectroMagneticField.h" + +namespace py = pybind11; + +namespace { +using DoubleArray = + py::array_t; + +std::vector to_vec(const DoubleArray &arr) { + auto buf = arr.request(); + return std::vector(static_cast(buf.ptr), + static_cast(buf.ptr) + buf.size); +} +} // namespace + +void init_GateMappedElectroMagneticField(py::module &m) { + + py::class_>( + m, "GateMappedElectroMagneticField") + + .def(py::init([](const G4VSolid *solid, + std::vector translations, + std::vector rotations, + double delta_chord_mm, int nx_B, int ny_B, int nz_B, + double x0_B, double y0_B, double z0_B, double dx_B, + double dy_B, double dz_B, DoubleArray Bx, DoubleArray By, + DoubleArray Bz, int nx_E, int ny_E, int nz_E, + double x0_E, double y0_E, double z0_E, double dx_E, + double dy_E, double dz_E, DoubleArray Ex, DoubleArray Ey, + DoubleArray Ez, + GateGridInterpolator::InterpolationMethod interp) { + GateGridInterpolator::GridDefinition gridDefB{ + nx_B, ny_B, nz_B, x0_B, y0_B, z0_B, dx_B, dy_B, dz_B}; + GateGridInterpolator::FieldValues fieldValuesB{ + to_vec(Bx), to_vec(By), to_vec(Bz)}; + GateGridInterpolator::GridDefinition gridDefE{ + nx_E, ny_E, nz_E, x0_E, y0_E, z0_E, dx_E, dy_E, dz_E}; + GateGridInterpolator::FieldValues fieldValuesE{ + to_vec(Ex), to_vec(Ey), to_vec(Ez)}; + return new GateMappedElectroMagneticField( + solid, translations, rotations, delta_chord_mm, gridDefB, + fieldValuesB, gridDefE, fieldValuesE, interp); + }), + py::arg("solid"), py::arg("translations"), py::arg("rotations"), + py::arg("delta_chord_mm"), py::arg("nx_B"), py::arg("ny_B"), + py::arg("nz_B"), py::arg("x0_B"), py::arg("y0_B"), py::arg("z0_B"), + py::arg("dx_B"), py::arg("dy_B"), py::arg("dz_B"), py::arg("Bx"), + py::arg("By"), py::arg("Bz"), py::arg("nx_E"), py::arg("ny_E"), + py::arg("nz_E"), py::arg("x0_E"), py::arg("y0_E"), py::arg("z0_E"), + py::arg("dx_E"), py::arg("dy_E"), py::arg("dz_E"), py::arg("Ex"), + py::arg("Ey"), py::arg("Ez"), py::arg("interpolation")) + + .def("SetTransforms", &GateMappedElectroMagneticField::SetTransforms, + py::arg("translations"), py::arg("rotations")); +} diff --git a/opengate/geometry/fields.py b/opengate/geometry/fields.py index 35c278153e..c977c341be 100644 --- a/opengate/geometry/fields.py +++ b/opengate/geometry/fields.py @@ -447,17 +447,17 @@ class UniformElectroMagneticField(ElectroMagneticField): """Uniform electromagnetic field with constant magnetic and electric field vectors.""" # hints for IDE - magnetic_field_vector: list - electric_field_vector: list + field_vector_B: list + field_vector_E: list user_info_defaults = { - "magnetic_field_vector": ( + "field_vector_B": ( [0, 0, 0], { "doc": "Magnetic field vector [Bx, By, Bz] in local volume coordinates.", }, ), - "electric_field_vector": ( + "field_vector_E": ( [0, 0, 0], { "doc": "Electric field vector [Ex, Ey, Ez] in local volume coordinates.", @@ -467,8 +467,8 @@ class UniformElectroMagneticField(ElectroMagneticField): def _create_inner_field(self): return g4.GateUniformElectroMagneticField( - g4.G4ThreeVector(*self.electric_field_vector), - g4.G4ThreeVector(*self.magnetic_field_vector), + g4.G4ThreeVector(*self.field_vector_E), + g4.G4ThreeVector(*self.field_vector_B), ) @@ -681,13 +681,105 @@ def create_field_manager(self, volume_obj) -> g4.G4FieldManager: class MappedElectroMagneticField(ElectroMagneticField): - """Electromagnetic field defined by values on a regular 3D Cartesian grid. + """Electromagnetic field with separate B and E grids on regular 3D Cartesian grids. - FIXME: implement this class. + field_matrix_B and field_matrix_E must have columns [x, y, z, Bx, By, Bz] + and [x, y, z, Ex, Ey, Ez] respectively, in Geant4 units. + Field lookup is performed entirely in C++ using trilinear or nearest-neighbour + interpolation, with B and E computed independently on their respective grids. """ + field_matrix_B: np.ndarray + field_matrix_E: np.ndarray + interpolation: str + + user_info_defaults = { + "field_matrix_B": ( + None, + { + "doc": ( + "2D numpy array on a regular Cartesian grid in Geant4 units. " + "Columns: [x, y, z, Bx, By, Bz]." + ), + }, + ), + "field_matrix_E": ( + None, + { + "doc": ( + "2D numpy array on a regular Cartesian grid in Geant4 units. " + "Columns: [x, y, z, Ex, Ey, Ez]." + ), + }, + ), + "interpolation": ( + "trilinear", + { + "doc": "Interpolation method: 'trilinear' (default) or 'nearest'.", + }, + ), + } + def create_field_manager(self, volume_obj) -> g4.G4FieldManager: - raise NotImplementedError("MappedElectroMagneticField is not yet implemented.") + if self.field_matrix_B is None: + raise ValueError( + "field_matrix_B must be provided for MappedElectroMagneticField" + ) + if self.field_matrix_E is None: + raise ValueError( + "field_matrix_E must be provided for MappedElectroMagneticField" + ) + if self.interpolation not in _interp_map: + raise ValueError( + f"Unknown interpolation '{self.interpolation}'. " + f"Choose 'trilinear' or 'nearest'." + ) + self._field_volume_obj = volume_obj + nx_B, ny_B, nz_B, x0_B, y0_B, z0_B, dx_B, dy_B, dz_B, (Bx, By, Bz) = ( + _parse_field_matrix( + self.field_matrix_B, "MappedElectroMagneticField (B grid)" + ) + ) + nx_E, ny_E, nz_E, x0_E, y0_E, z0_E, dx_E, dy_E, dz_E, (Ex, Ey, Ez) = ( + _parse_field_matrix( + self.field_matrix_E, "MappedElectroMagneticField (E grid)" + ) + ) + g4_translations, g4_rotations = self._make_g4_transforms() + gate_field = g4.GateMappedElectroMagneticField( + self._field_volume_obj.g4_solid, + g4_translations, + g4_rotations, + self.delta_chord, + nx_B, + ny_B, + nz_B, + x0_B, + y0_B, + z0_B, + dx_B, + dy_B, + dz_B, + Bx, + By, + Bz, + nx_E, + ny_E, + nz_E, + x0_E, + y0_E, + z0_E, + dx_E, + dy_E, + dz_E, + Ex, + Ey, + Ez, + _interp_map[self.interpolation], + ) + return self._build_field_manager( + None, gate_field, g4.G4EqMagElectricField, 8, volume_obj + ) field_types = { From e8b4ef98cb32078b61bbecff47ff14106c082515 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Thu, 7 May 2026 08:02:01 +0200 Subject: [PATCH 15/28] Refactor tests and add new ones --- opengate/geometry/fields.py | 2 +- .../geometry/test099_fields_analytical_B.py | 2 +- .../geometry/test099_fields_analytical_E.py | 2 +- .../tests/src/geometry/test099_fields_api.py | 10 +- .../test099_fields_custom_vs_native_B.py | 12 +- .../test099_fields_custom_vs_native_E.py | 12 +- ...t099_fields_mapped_multi_volume_refresh.py | 177 ++++++++++++++++ ...st099_fields_mapped_repeated_placements.py | 198 ++++++++++++++++++ .../test099_fields_mapped_rotated_volume.py | 156 ++++++++++++++ .../test099_fields_mapped_vs_uniform_B.py | 10 +- .../test099_fields_mapped_vs_uniform_E.py | 2 +- .../geometry/test099_fields_serialization.py | 2 +- 12 files changed, 558 insertions(+), 27 deletions(-) create mode 100644 opengate/tests/src/geometry/test099_fields_mapped_multi_volume_refresh.py create mode 100644 opengate/tests/src/geometry/test099_fields_mapped_repeated_placements.py create mode 100644 opengate/tests/src/geometry/test099_fields_mapped_rotated_volume.py diff --git a/opengate/geometry/fields.py b/opengate/geometry/fields.py index c977c341be..04cd6821b9 100644 --- a/opengate/geometry/fields.py +++ b/opengate/geometry/fields.py @@ -84,7 +84,7 @@ def __init__(self, *args, **kwargs) -> None: self.attached_to = [] self._field_changes_energy = False - # Integration objects — shared by all field subclasses + # Integration objects - shared by all field subclasses self.g4_equation_of_motion = None self.g4_integrator_stepper = None self.g4_integration_driver = None diff --git a/opengate/tests/src/geometry/test099_fields_analytical_B.py b/opengate/tests/src/geometry/test099_fields_analytical_B.py index 39415c484f..c78207fb65 100644 --- a/opengate/tests/src/geometry/test099_fields_analytical_B.py +++ b/opengate/tests/src/geometry/test099_fields_analytical_B.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 """ -Test 099 — Uniform magnetic field: analytical validation. +Test 099 - Uniform magnetic field: analytical validation. A proton enters a region with uniform B along Y. Checks: diff --git a/opengate/tests/src/geometry/test099_fields_analytical_E.py b/opengate/tests/src/geometry/test099_fields_analytical_E.py index 9203b5d286..f97ab4927b 100644 --- a/opengate/tests/src/geometry/test099_fields_analytical_E.py +++ b/opengate/tests/src/geometry/test099_fields_analytical_E.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 """ -Test 099 — Uniform electric field: analytical validation. +Test 099 - Uniform electric field: analytical validation. A proton enters a region with uniform E along X. Checks: diff --git a/opengate/tests/src/geometry/test099_fields_api.py b/opengate/tests/src/geometry/test099_fields_api.py index 95bd3f1316..c5e1df4508 100644 --- a/opengate/tests/src/geometry/test099_fields_api.py +++ b/opengate/tests/src/geometry/test099_fields_api.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 """ -Test 099 — Field API tests. +Test 099 - Field API tests. Checks the following behaviours: 1. add_field() before volume is added to simulation -> fatal @@ -23,10 +23,10 @@ def _expect_error(fn, error_msg): """Run fn, expect an exception (fatal() raises Exception, others raise specific types).""" try: fn() - print(f" FAIL — expected error: {error_msg}") + print(f" FAIL - expected error: {error_msg}") return False except Exception as e: - print(f" OK — got {type(e).__name__}: {error_msg}") + print(f" OK - got {type(e).__name__}: {error_msg}") return True @@ -120,9 +120,9 @@ def _expect_error(fn, error_msg): and box1.field == "B_shared" and box2.field == "B_shared" ) - print(f" {'OK' if ok else 'FAIL'} — same field on two volumes") + print(f" {'OK' if ok else 'FAIL'} - same field on two volumes") except Exception as e: - print(f" FAIL — unexpected error: {e}") + print(f" FAIL - unexpected error: {e}") ok = False is_ok = is_ok and ok diff --git a/opengate/tests/src/geometry/test099_fields_custom_vs_native_B.py b/opengate/tests/src/geometry/test099_fields_custom_vs_native_B.py index d2db37b871..0edc9ccf43 100644 --- a/opengate/tests/src/geometry/test099_fields_custom_vs_native_B.py +++ b/opengate/tests/src/geometry/test099_fields_custom_vs_native_B.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 """ -Test 099 — Custom trampoline B field vs native G4 uniform B field. +Test 099 - Custom trampoline B field vs native G4 uniform B field. Places two boxes side-by-side in one simulation: one with the native G4UniformMagField, the other with a CustomMagneticField returning the @@ -141,10 +141,10 @@ def custom_uniform_B(x, y, z, t): is_ok_z = np.all(dz < r_TOL) is_ok_e = np.all(dKE < e_TOL) - print(f"dx max: {dx.max():.6f} mm — OK: {is_ok_x}") - print(f"dy max: {dy.max():.6f} mm — OK: {is_ok_y}") - print(f"dz max: {dz.max():.6f} mm — OK: {is_ok_z}") - print(f"dKE max: {dKE.max():.6f} MeV — OK: {is_ok_e}") + print(f"dx max: {dx.max():.6f} mm - OK: {is_ok_x}") + print(f"dy max: {dy.max():.6f} mm - OK: {is_ok_y}") + print(f"dz max: {dz.max():.6f} mm - OK: {is_ok_z}") + print(f"dKE max: {dKE.max():.6f} MeV - OK: {is_ok_e}") # Compare custom to analytical r = cyclotron_radius(T, By, PROTON_MASS, 1 * g4_eplus) @@ -152,7 +152,7 @@ def custom_uniform_B(x, y, z, t): cx, cz = 0.0 - r, -box_half_z residual = np.sqrt((custom_x - cx) ** 2 + (custom_z - cz) ** 2) - r is_ok_analytical = np.all(np.abs(residual) < r_TOL) - print(f"Custom vs analytical circle residual: {residual} — OK: {is_ok_analytical}") + print(f"Custom vs analytical circle residual: {residual} - OK: {is_ok_analytical}") is_ok = is_ok_x and is_ok_y and is_ok_z and is_ok_e and is_ok_analytical utility.test_ok(is_ok) diff --git a/opengate/tests/src/geometry/test099_fields_custom_vs_native_E.py b/opengate/tests/src/geometry/test099_fields_custom_vs_native_E.py index 0492b049cc..09042ebeee 100644 --- a/opengate/tests/src/geometry/test099_fields_custom_vs_native_E.py +++ b/opengate/tests/src/geometry/test099_fields_custom_vs_native_E.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 """ -Test 099 — Custom trampoline E field vs native G4 uniform E field. +Test 099 - Custom trampoline E field vs native G4 uniform E field. Places two boxes side-by-side in one simulation: one with the native G4UniformElectricField, the other with a CustomElectricField returning the @@ -140,16 +140,16 @@ def custom_uniform_E(x, y, z, t): is_ok_z = np.all(dz < r_TOL) is_ok_e = np.all(dKE < e_TOL) - print(f"dx max: {dx.max():.6f} mm — OK: {is_ok_x}") - print(f"dy max: {dy.max():.6f} mm — OK: {is_ok_y}") - print(f"dz max: {dz.max():.6f} mm — OK: {is_ok_z}") - print(f"dKE max: {dKE.max():.6f} MeV — OK: {is_ok_e}") + print(f"dx max: {dx.max():.6f} mm - OK: {is_ok_x}") + print(f"dy max: {dy.max():.6f} mm - OK: {is_ok_y}") + print(f"dz max: {dz.max():.6f} mm - OK: {is_ok_z}") + print(f"dKE max: {dKE.max():.6f} MeV - OK: {is_ok_e}") # Compare custom to analytical x_from_energy = (custom_KE - T) / (q * Ex) is_ok_analytical = np.all(np.abs(custom_x - x_from_energy) < r_TOL) print( - f"Custom vs analytical x residual: {custom_x - x_from_energy} — OK: {is_ok_analytical}" + f"Custom vs analytical x residual: {custom_x - x_from_energy} - OK: {is_ok_analytical}" ) is_ok = is_ok_x and is_ok_y and is_ok_z and is_ok_e and is_ok_analytical diff --git a/opengate/tests/src/geometry/test099_fields_mapped_multi_volume_refresh.py b/opengate/tests/src/geometry/test099_fields_mapped_multi_volume_refresh.py new file mode 100644 index 0000000000..897f4f8eaf --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_mapped_multi_volume_refresh.py @@ -0,0 +1,177 @@ +#!/usr/bin/env python3 +""" +Test 099 - MappedMagneticField attached to two volumes, dynamic geometry between runs. + +Same scenario as test099_fields_multi_volume_refresh, but the shared field is a +MappedMagneticField. +""" + +import itertools + +import numpy as np +import uproot +from scipy.spatial.transform import Rotation + +import opengate as gate +from opengate.geometry import fields +from opengate.tests import utility +from test099_fields_helpers import ( + g4_mm, + g4_MeV, + g4_tesla, + g4_eplus, + PROTON_MASS, + magnetic_deflection, +) + +g4_s = gate.g4_units.s + +BOX_SIZE = 100 * g4_mm +SEP = 200 * g4_mm +By = 3 * g4_tesla +KE = 10 * g4_MeV + +EXPECTED_DEFLECTION = magnetic_deflection(KE, By, PROTON_MASS, 1 * g4_eplus, BOX_SIZE) +TOL_DEFLECTED = 0.01 * g4_mm +TOL_ZERO = 0.01 * g4_mm + +VISU = False + + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, output_folder="test099_fields") + + sim = gate.Simulation() + sim.g4_verbose = False + sim.visu = False + sim.random_seed = 99001 + sim.number_of_threads = 1 + sim.output_dir = paths.output + + if VISU: + sim.visu = True + sim.visu_type = "qt" + sim.visu_commands.append("/vis/scene/endOfEventAction accumulate") + sim.visu_commands.append("/vis/scene/add/trajectories smooth") + sim.visu_commands.append("/vis/scene/add/magneticField 20 fullArrow") + + sim.run_timing_intervals = [ + [0, 0.5 * g4_s], + [0.5 * g4_s, 1.0 * g4_s], + ] + + sim.world.size = [1 * gate.g4_units.m, 1 * gate.g4_units.m, 1 * gate.g4_units.m] + sim.world.material = "G4_Galactic" + + box1 = sim.add_volume("Box", "box1") + box1.size = [BOX_SIZE, BOX_SIZE, BOX_SIZE] + box1.material = "G4_Galactic" + box1.translation = [-SEP, 0, 0] + box1.add_dynamic_parametrisation( + rotation=[ + Rotation.from_euler("z", 0, degrees=True).as_matrix(), + Rotation.from_euler("z", 90, degrees=True).as_matrix(), + ] + ) + + box2 = sim.add_volume("Box", "box2") + box2.size = [BOX_SIZE, BOX_SIZE, BOX_SIZE] + box2.material = "G4_Galactic" + box2.translation = [+SEP, 0, 0] + + half = BOX_SIZE / 2 + corners = list(itertools.product([-half, half], repeat=3)) + field_matrix = np.array([[x, y, z, 0.0, By, 0.0] for x, y, z in corners]) + + field = fields.MappedMagneticField(name="B_mapped") + field.field_matrix = field_matrix + field.interpolation = "trilinear" + box1.add_field(field) + box2.add_field(field) + + src1 = sim.add_source("GenericSource", "src_box1") + src1.particle = "proton" + src1.n = [1, 1] + src1.energy.type = "mono" + src1.energy.mono = KE + src1.position.type = "point" + src1.position.translation = [-SEP, 0, -300 * g4_mm] + src1.direction.type = "momentum" + src1.direction.momentum = [0, 0, 1] + + src2 = sim.add_source("GenericSource", "src_box2") + src2.particle = "proton" + src2.n = [1, 1] + src2.energy.type = "mono" + src2.energy.mono = KE + src2.position.type = "point" + src2.position.translation = [+SEP, 0, -300 * g4_mm] + src2.direction.type = "momentum" + src2.direction.momentum = [0, 0, 1] + + phsp1 = sim.add_actor("PhaseSpaceActor", "phsp_box1") + phsp1.attached_to = "box1" + phsp1.attributes = ["PostPosition", "RunID"] + phsp1.output_filename = "phsp_mapped_refresh_box1.root" + phsp1.steps_to_store = "exiting" + phsp1.root_output.write_to_disk = True + + phsp2 = sim.add_actor("PhaseSpaceActor", "phsp_box2") + phsp2.attached_to = "box2" + phsp2.attributes = ["PostPosition", "RunID"] + phsp2.output_filename = "phsp_mapped_refresh_box2.root" + phsp2.steps_to_store = "exiting" + phsp2.root_output.write_to_disk = True + + sim.run() + + df1 = uproot.open(str(paths.output / "phsp_mapped_refresh_box1.root"))[ + "phsp_box1;1" + ].arrays(library="pd") + df2 = uproot.open(str(paths.output / "phsp_mapped_refresh_box2.root"))[ + "phsp_box2;1" + ].arrays(library="pd") + + r0b1 = df1[df1["RunID"] == 0].iloc[0] + r1b1 = df1[df1["RunID"] == 1].iloc[0] + r0b2 = df2[df2["RunID"] == 0].iloc[0] + r1b2 = df2[df2["RunID"] == 1].iloc[0] + + x1_r0 = r0b1["PostPosition_X"] - (-SEP) + x1_r1 = r1b1["PostPosition_X"] - (-SEP) + y1_r0 = r0b1["PostPosition_Y"] + y1_r1 = r1b1["PostPosition_Y"] + + x2_r0 = r0b2["PostPosition_X"] - (+SEP) + x2_r1 = r1b2["PostPosition_X"] - (+SEP) + y2_r0 = r0b2["PostPosition_Y"] + y2_r1 = r1b2["PostPosition_Y"] + + print( + f"Analytical expected deflection: {EXPECTED_DEFLECTION:.2f} mm (threshold: {TOL_DEFLECTED:.2f} mm)" + ) + print(f"box1 run0: ΔX={x1_r0:.2f} Y={y1_r0:.2f} mm") + print(f"box1 run1: ΔX={x1_r1:.2f} Y={y1_r1:.2f} mm") + print(f"box2 run0: ΔX={x2_r0:.2f} Y={y2_r0:.2f} mm") + print(f"box2 run1: ΔX={x2_r1:.2f} Y={y2_r1:.2f} mm") + + ok_r0b1 = ( + np.abs(x1_r0 + EXPECTED_DEFLECTION) < TOL_DEFLECTED and abs(y1_r0) < TOL_ZERO + ) + ok_r0b2 = ( + np.abs(x2_r0 + EXPECTED_DEFLECTION) < TOL_DEFLECTED and abs(y2_r0) < TOL_ZERO + ) + ok_r1b2 = ( + np.abs(x2_r1 + EXPECTED_DEFLECTION) < TOL_DEFLECTED and abs(y2_r1) < TOL_ZERO + ) + ok_r1b1 = ( + np.abs(y1_r1 + EXPECTED_DEFLECTION) < TOL_DEFLECTED and abs(x1_r1) < TOL_ZERO + ) + + print(f"\nbox1 run0 deflects in -X: {ok_r0b1}") + print(f"box2 run0 deflects in -X: {ok_r0b2}") + print(f"box2 run1 deflects in -X (unchanged): {ok_r1b2}") + print(f"box1 run1 deflects in -Y (verifies refresh): {ok_r1b1}") + + is_ok = ok_r0b1 and ok_r0b2 and ok_r1b2 and ok_r1b1 + utility.test_ok(is_ok) diff --git a/opengate/tests/src/geometry/test099_fields_mapped_repeated_placements.py b/opengate/tests/src/geometry/test099_fields_mapped_repeated_placements.py new file mode 100644 index 0000000000..f67cef4bf9 --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_mapped_repeated_placements.py @@ -0,0 +1,198 @@ +#!/usr/bin/env python3 +""" +Test 099 - MappedMagneticField on a repeated volume. + +Same scenario as test099_fields_repeated_placements, but both fields are +MappedMagneticField. +""" + +import itertools + +import numpy as np +import uproot + +import opengate as gate +from opengate.geometry import fields +from opengate.geometry.utility import get_grid_repetition +from opengate.tests import utility +from test099_fields_helpers import ( + g4_mm, + g4_MeV, + g4_tesla, + g4_eplus, + PROTON_MASS, + magnetic_deflection, +) + +BOX_SIZE = 50 * g4_mm +N_REPS = 10 +SLAB_SIZE = BOX_SIZE / N_REPS +SPACING = SLAB_SIZE +By = 3 * g4_tesla +KE = 10 * g4_MeV +SEP = 100 * g4_mm + +EXPECTED_DEFLECTION = magnetic_deflection(KE, By, PROTON_MASS, 1 * g4_eplus, BOX_SIZE) + +TOL_DEFLECTED = 0.1 * g4_mm +TOL_Y = 0.1 * g4_mm +TOL_E = 0.01 * g4_MeV + +VISU = False + + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, output_folder="test099_fields") + + sim = gate.Simulation() + sim.g4_verbose = False + sim.visu = False + sim.random_seed = 12345 + sim.number_of_threads = 1 + sim.output_dir = paths.output + + if VISU: + sim.visu = True + sim.visu_type = "qt" + sim.visu_commands.append("/vis/scene/endOfEventAction accumulate") + sim.visu_commands.append("/vis/scene/add/trajectories smooth") + sim.visu_commands.append("/vis/scene/add/magneticField 50 fullArrow") + + m = gate.g4_units.m + sim.world.size = [2 * m, 1 * m, 1 * m] + sim.world.material = "G4_Galactic" + + # --- LEFT: single box at X = -SEP --- + box_single = sim.add_volume("Box", "box_single") + box_single.size = [BOX_SIZE, BOX_SIZE, BOX_SIZE] + box_single.material = "G4_Galactic" + box_single.translation = [-SEP, 0, 0] + + # 2x2x2 grid covering the full single-box local extent. + half = BOX_SIZE / 2 + corners_single = list(itertools.product([-half, half], repeat=3)) + mat_single = np.array([[x, y, z, 0.0, By, 0.0] for x, y, z in corners_single]) + + field_single = fields.MappedMagneticField(name="B_single") + field_single.field_matrix = mat_single + field_single.interpolation = "trilinear" + box_single.add_field(field_single) + + src_single = sim.add_source("GenericSource", "src_single") + src_single.particle = "proton" + src_single.n = 1 + src_single.energy.type = "mono" + src_single.energy.mono = KE + src_single.position.type = "point" + src_single.position.translation = [-SEP, 0, -200 * g4_mm] + src_single.direction.type = "momentum" + src_single.direction.momentum = [0, 0, 1] + + phsp_single = sim.add_actor("PhaseSpaceActor", "phsp_single") + phsp_single.attached_to = "box_single" + phsp_single.attributes = ["PostKineticEnergy", "PostPosition"] + phsp_single.output_filename = "phsp_mapped_repeated_single.root" + phsp_single.steps_to_store = "exiting" + phsp_single.root_output.write_to_disk = True + + # --- RIGHT: repeated box (N_REPS slabs along Z) at X = +SEP --- + translations = get_grid_repetition( + [1, 1, N_REPS], [0, 0, SPACING], start=[0, 0, -SPACING] + ) + translations = [[+SEP + t[0], t[1], t[2]] for t in translations] + + box_multi = sim.add_volume("Box", "box_multi") + box_multi.size = [BOX_SIZE, BOX_SIZE, SLAB_SIZE] + box_multi.material = "G4_Galactic" + box_multi.translation = translations + + # 2x2x2 grid covering the slab local extent: BOX_SIZE x BOX_SIZE x SLAB_SIZE. + hxy = BOX_SIZE / 2 + hz = SLAB_SIZE / 2 + corners_slab = list(itertools.product([-hxy, hxy], [-hxy, hxy], [-hz, hz])) + mat_multi = np.array([[x, y, z, 0.0, By, 0.0] for x, y, z in corners_slab]) + + field_multi = fields.MappedMagneticField(name="B_multi") + field_multi.field_matrix = mat_multi + field_multi.interpolation = "trilinear" + box_multi.add_field(field_multi) + + src_multi = sim.add_source("GenericSource", "src_multi") + src_multi.particle = "proton" + src_multi.n = 1 + src_multi.energy.type = "mono" + src_multi.energy.mono = KE + src_multi.position.type = "point" + src_multi.position.translation = [+SEP, 0, -200 * g4_mm] + src_multi.direction.type = "momentum" + src_multi.direction.momentum = [0, 0, 1] + + phsp_multi = sim.add_actor("PhaseSpaceActor", "phsp_multi") + phsp_multi.attached_to = "box_multi" + phsp_multi.attributes = ["PostKineticEnergy", "PostPosition"] + phsp_multi.output_filename = "phsp_mapped_repeated_multi.root" + phsp_multi.steps_to_store = "exiting" + phsp_multi.root_output.write_to_disk = True + + sim.run() + + df_single = uproot.open(str(paths.output / "phsp_mapped_repeated_single.root"))[ + "phsp_single;1" + ].arrays(library="pd") + df_multi = uproot.open(str(paths.output / "phsp_mapped_repeated_multi.root"))[ + "phsp_multi;1" + ].arrays(library="pd") + + row_single = df_single.sort_values("PostPosition_Z").iloc[-1] + row_multi = df_multi.sort_values("PostPosition_Z").iloc[-1] + + x_single = row_single["PostPosition_X"] - (-SEP) + y_single = row_single["PostPosition_Y"] + ke_single = row_single["PostKineticEnergy"] + + x_multi = row_multi["PostPosition_X"] - (+SEP) + y_multi = row_multi["PostPosition_Y"] + ke_multi = row_multi["PostKineticEnergy"] + + print( + f"Analytical expected deflection: {-EXPECTED_DEFLECTION:.2f} mm (tolerance: ±{TOL_DEFLECTED:.2f} mm)" + ) + print( + f"Single box ({BOX_SIZE:.0f} mm): ΔX={x_single:.2f} mm Y={y_single:.2f} mm KE={ke_single:.4f} MeV" + ) + print( + f"{N_REPS}x slabs ({SLAB_SIZE:.1f} mm each): ΔX={x_multi:.2f} mm Y={y_multi:.2f} mm KE={ke_multi:.4f} MeV" + ) + + ok_single = abs(x_single + EXPECTED_DEFLECTION) < TOL_DEFLECTED + ok_multi = abs(x_multi + EXPECTED_DEFLECTION) < TOL_DEFLECTED + print( + f"Single box matches analytical: {ok_single} (|{x_single:.2f} - {-EXPECTED_DEFLECTION:.2f}| < {TOL_DEFLECTED:.2f})" + ) + print( + f"{N_REPS}-slabs match analytical: {ok_multi} (|{x_multi:.2f} - {-EXPECTED_DEFLECTION:.2f}| < {TOL_DEFLECTED:.2f})" + ) + + ok_y_single = abs(y_single) < TOL_Y + ok_y_multi = abs(y_multi) < TOL_Y + print(f"No Y deflection (single): {ok_y_single} ({y_single:.3f} mm)") + print(f"No Y deflection (multi): {ok_y_multi} ({y_multi:.3f} mm)") + + ok_energy_single = abs(ke_single - KE) < TOL_E + ok_energy_multi = abs(ke_multi - KE) < TOL_E + print( + f"Energy conserved (single): {ok_energy_single} (dKE={ke_single - KE:.4f} MeV)" + ) + print( + f"Energy conserved (multi): {ok_energy_multi} (dKE={ke_multi - KE:.4f} MeV)" + ) + + is_ok = ( + ok_single + and ok_multi + and ok_y_single + and ok_y_multi + and ok_energy_single + and ok_energy_multi + ) + utility.test_ok(is_ok) diff --git a/opengate/tests/src/geometry/test099_fields_mapped_rotated_volume.py b/opengate/tests/src/geometry/test099_fields_mapped_rotated_volume.py new file mode 100644 index 0000000000..00433529b0 --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_mapped_rotated_volume.py @@ -0,0 +1,156 @@ +#!/usr/bin/env python3 +""" +Test 099 - MappedMagneticField direction follows volume rotation. + +Same scenario as test099_fields_rotated_volume, but the shared field is a +MappedMagneticField. +""" + +import itertools + +import numpy as np +import uproot +from scipy.spatial.transform import Rotation + +import opengate as gate +from opengate.geometry import fields +from opengate.tests import utility +from test099_fields_helpers import ( + g4_mm, + g4_MeV, + g4_tesla, + g4_eplus, + PROTON_MASS, + magnetic_deflection, +) + +BOX_SIZE = 100 * g4_mm +By = 3 * g4_tesla +KE = 10 * g4_MeV +SEP = 300 * g4_mm + +EXPECTED_DEFLECTION = magnetic_deflection(KE, By, PROTON_MASS, 1 * g4_eplus, BOX_SIZE) +TOL_DEFLECTED = 0.01 * g4_mm +TOL_ZERO = 0.01 * g4_mm + +VISU = False + + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, output_folder="test099_fields") + + sim = gate.Simulation() + sim.g4_verbose = False + sim.visu = False + sim.random_seed = 42 + sim.number_of_threads = 1 + sim.output_dir = paths.output + + if VISU: + sim.visu = True + sim.visu_type = "qt" + sim.visu_commands.append("/vis/scene/endOfEventAction accumulate") + sim.visu_commands.append("/vis/scene/add/trajectories smooth") + sim.visu_commands.append("/vis/scene/add/magneticField 20 fullArrow") + + m = gate.g4_units.m + sim.world.size = [2 * m, 2 * m, 2 * m] + sim.world.material = "G4_Galactic" + + box_a = sim.add_volume("Box", "box_a") + box_a.size = [BOX_SIZE, BOX_SIZE, BOX_SIZE] + box_a.material = "G4_Galactic" + box_a.translation = [-SEP, 0, 0] + + box_b = sim.add_volume("Box", "box_b") + box_b.size = [BOX_SIZE, BOX_SIZE, BOX_SIZE] + box_b.material = "G4_Galactic" + box_b.translation = [+SEP, 0, 0] + box_b.rotation = Rotation.from_euler("z", 90, degrees=True).as_matrix() + + half = BOX_SIZE / 2 + corners = list(itertools.product([-half, half], repeat=3)) + field_matrix = np.array([[x, y, z, 0.0, By, 0.0] for x, y, z in corners]) + + field = fields.MappedMagneticField(name="B_mapped") + field.field_matrix = field_matrix + field.interpolation = "trilinear" + box_a.add_field(field) + box_b.add_field(field) + + src_a = sim.add_source("GenericSource", "src_a") + src_a.particle = "proton" + src_a.n = 1 + src_a.energy.type = "mono" + src_a.energy.mono = KE + src_a.position.type = "point" + src_a.position.translation = [-SEP, 0, -300 * g4_mm] + src_a.direction.type = "momentum" + src_a.direction.momentum = [0, 0, 1] + + src_b = sim.add_source("GenericSource", "src_b") + src_b.particle = "proton" + src_b.n = 1 + src_b.energy.type = "mono" + src_b.energy.mono = KE + src_b.position.type = "point" + src_b.position.translation = [+SEP, 0, -300 * g4_mm] + src_b.direction.type = "momentum" + src_b.direction.momentum = [0, 0, 1] + + phsp_a = sim.add_actor("PhaseSpaceActor", "phsp_a") + phsp_a.attached_to = "box_a" + phsp_a.attributes = ["PostPosition"] + phsp_a.output_filename = "phsp_mapped_rotated_a.root" + phsp_a.steps_to_store = "exiting" + phsp_a.root_output.write_to_disk = True + + phsp_b = sim.add_actor("PhaseSpaceActor", "phsp_b") + phsp_b.attached_to = "box_b" + phsp_b.attributes = ["PostPosition"] + phsp_b.output_filename = "phsp_mapped_rotated_b.root" + phsp_b.steps_to_store = "exiting" + phsp_b.root_output.write_to_disk = True + + sim.run() + + df_a = uproot.open(str(paths.output / "phsp_mapped_rotated_a.root"))[ + "phsp_a;1" + ].arrays(library="pd") + df_b = uproot.open(str(paths.output / "phsp_mapped_rotated_b.root"))[ + "phsp_b;1" + ].arrays(library="pd") + + row_a = df_a.sort_values("PostPosition_Z").iloc[-1] + row_b = df_b.sort_values("PostPosition_Z").iloc[-1] + + xa = row_a["PostPosition_X"] - (-SEP) + ya = row_a["PostPosition_Y"] + xb = row_b["PostPosition_X"] - (+SEP) + yb = row_b["PostPosition_Y"] + + print( + f"Analytical expected deflection: {-EXPECTED_DEFLECTION:.2f} mm (threshold: {TOL_DEFLECTED:.2f} mm)" + ) + print(f"box_a (unrotated): X={xa:.2f} mm Y={ya:.2f} mm") + print(f"box_b (R_z 90°): X={xb:.2f} mm Y={yb:.2f} mm") + + ok_a_x = np.abs(xa + EXPECTED_DEFLECTION) < TOL_DEFLECTED + ok_a_y = np.abs(ya) < TOL_ZERO + print( + f"box_a deflects in -X: {ok_a_x} (|{xa:.2f} - {-EXPECTED_DEFLECTION:.2f}| < {TOL_DEFLECTED:.2f})" + ) + print(f"box_a no Y deflection: {ok_a_y} ({ya:.2f} mm)") + + ok_b_y = np.abs(yb + EXPECTED_DEFLECTION) < TOL_DEFLECTED + ok_b_x = np.abs(xb) < TOL_ZERO + print( + f"box_b deflects in -Y: {ok_b_y} (|{yb:.2f} - {-EXPECTED_DEFLECTION:.2f}| < {TOL_DEFLECTED:.2f})" + ) + print(f"box_b no X deflection: {ok_b_x} ({xb:.2f} mm)") + + ok_symmetry = np.abs(xa - yb) < TOL_DEFLECTED + print(f"box_a deflection in -X is similar to box_b deflection in -Y: {ok_symmetry}") + + is_ok = ok_a_x and ok_a_y and ok_b_y and ok_b_x and ok_symmetry + utility.test_ok(is_ok) diff --git a/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_B.py b/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_B.py index a98668051f..c094cd8b7f 100644 --- a/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_B.py +++ b/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_B.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 """ -Test 099 — MappedMagneticField with uniform grid vs UniformMagneticField. +Test 099 - MappedMagneticField with uniform grid vs UniformMagneticField. Places two boxes side-by-side: one with UniformMagneticField, the other with a MappedMagneticField whose grid contains a constant B value everywhere. A @@ -141,9 +141,9 @@ is_ok_z = np.all(dz < r_TOL) is_ok_e = np.all(dKE < e_TOL) - print(f"Uniform vs mapped — dx max : {dx.max() / g4_mm:.6f} mm — OK: {is_ok_x}") - print(f"Uniform vs mapped — dz max : {dz.max() / g4_mm:.6f} mm — OK: {is_ok_z}") - print(f"Uniform vs mapped — dKE max: {dKE.max() / g4_MeV:.6f} MeV — OK: {is_ok_e}") + print(f"Uniform vs mapped - dx max : {dx.max() / g4_mm:.6f} mm - OK: {is_ok_x}") + print(f"Uniform vs mapped - dz max : {dz.max() / g4_mm:.6f} mm - OK: {is_ok_z}") + print(f"Uniform vs mapped - dKE max: {dKE.max() / g4_MeV:.6f} MeV - OK: {is_ok_e}") # Analytical cyclotron check on the mapped result r = cyclotron_radius(T, By, PROTON_MASS, 1 * g4_eplus) @@ -151,7 +151,7 @@ residual = np.sqrt((mapped_x - cx) ** 2 + (mapped_z - cz) ** 2) - r is_ok_analytical = np.all(np.abs(residual) < r_TOL) print( - f"Mapped vs analytical circle residual: {residual / g4_mm} mm — OK: {is_ok_analytical}" + f"Mapped vs analytical circle residual: {residual / g4_mm} mm - OK: {is_ok_analytical}" ) is_ok = is_ok_x and is_ok_z and is_ok_e and is_ok_analytical diff --git a/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py b/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py index a0129bf815..3931a08e42 100644 --- a/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py +++ b/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 """ -Test 099 — MappedMagneticField with uniform grid vs UniformMagneticField. +Test 099 - MappedMagneticField with uniform grid vs UniformMagneticField. Places two boxes side-by-side: one with UniformMagneticField, the other with a MappedMagneticField whose grid contains a constant B value everywhere. A diff --git a/opengate/tests/src/geometry/test099_fields_serialization.py b/opengate/tests/src/geometry/test099_fields_serialization.py index d3c80adb5b..531f3e7f22 100644 --- a/opengate/tests/src/geometry/test099_fields_serialization.py +++ b/opengate/tests/src/geometry/test099_fields_serialization.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 """ -Test 099 — Field serialization round-trip. +Test 099 - Field serialization round-trip. Tests that non-custom field types survive serialization and that custom fields correctly refuse it. From 7a00fdb0011e8600cd441c96ccda299f299f749c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Thu, 7 May 2026 08:17:36 +0200 Subject: [PATCH 16/28] Documentation for mapped fields + fix related build warnings --- docs/source/conf.py | 5 ++ docs/source/user_guide/user_guide_fields.rst | 78 ++++++++++++++++++-- 2 files changed, 78 insertions(+), 5 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index 9b004847b4..c0457cc407 100755 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -128,6 +128,11 @@ # Enable nitpicky mode for broken references/links nitpicky = True +# Suppress warnings for pybind11-bound C++ types that have no Sphinx docs entry +nitpick_ignore = [ + ("py:class", "opengate_core.G4FieldManager"), +] + # -- Options for HTML output ------------------------------------------------- # The theme to use for HTML and HTML Help pages. See the documentation for diff --git a/docs/source/user_guide/user_guide_fields.rst b/docs/source/user_guide/user_guide_fields.rst index b8cfa8e2d6..5081909b55 100644 --- a/docs/source/user_guide/user_guide_fields.rst +++ b/docs/source/user_guide/user_guide_fields.rst @@ -76,6 +76,30 @@ Magnetic fields Custom fields are evaluated via a Python callback for every field evaluation during tracking. This will significantly slow down the simulation. Prefer native types when possible. We are currently working on developing a faster custom field implementation. For more details, see :ref:`user_guide_fields_performance`. +**MappedMagneticField** -- Magnetic field defined by values on a regular 3D Cartesian grid. Field values are interpolated between grid points (trilinear by default). This is the recommended approach for importing fields from external calculations (e.g. finite element solvers) or for replacing a slow ``CustomMagneticField`` with a pre-sampled C++ equivalent. + +The field is specified as a 2D array with columns ``[x, y, z, Bx, By, Bz]`` in Geant4 internal units. The grid must be regular (uniform spacing along each axis) and complete (every combination of the sampled x, y, z values must be present). All coordinates are in the local frame of the attached volume. Degenerate axes are not allowed, so the minimum valid grid is 2×2×2 (eight corner points). + +.. code-block:: python + + import numpy as np + + # data.csv has columns: x, y, z (in mm), Bx, By, Bz (in T) + mm = gate.g4_units.mm + tesla = gate.g4_units.tesla + + field_matrix = np.loadtxt("data.csv", delimiter=",") + field_matrix[:, :3] *= mm # convert positions to Geant4 internal length units + field_matrix[:, 3:] *= tesla # convert field to Geant4 internal field units + + field = fields.MappedMagneticField(name="B_mapped") + field.field_matrix = field_matrix + box.add_field(field) + +.. warning:: Points outside the grid + + If the grid does not cover the entire volume, field values will be extrapolated outside the grid using the nearest valid value (i.e. clamping to the edge). We recommend defining the grid to cover the entire volume to avoid unexpected behaviour. + Electric fields ~~~~~~~~~~~~~~~ @@ -102,6 +126,18 @@ Electric fields field.field_function = my_E_field box.add_field(field) +.. warning:: Performance warning + + Same GIL overhead as ``CustomMagneticField``. See :ref:`user_guide_fields_performance`. + +**MappedElectricField** -- Electric field defined on a regular 3D Cartesian grid, same interface as ``MappedMagneticField``. Columns are ``[x, y, z, Ex, Ey, Ez]``. + +.. code-block:: python + + field = fields.MappedElectricField(name="E_mapped") + field.field_matrix = field_matrix # columns: [x, y, z, Ex, Ey, Ez], values in Geant4 internal units + box.add_field(field) + Electromagnetic fields ~~~~~~~~~~~~~~~~~~~~~~ @@ -116,7 +152,8 @@ Combined magnetic and electric fields. m = gate.g4_units.m field = fields.UniformElectroMagneticField(name="EM_uniform") - field.field_vector = [0, 0, 1 * tesla, 1e6 * volt / m, 0, 0] # [Bx, By, Bz, Ex, Ey, Ez] <- relative to the volume's local coordinate system + field.field_vector_B = [0, 0, 1 * tesla] # [Bx, By, Bz] in local coordinates + field.field_vector_E = [1e6 * volt / m, 0, 0] # [Ex, Ey, Ez] in local coordinates box.add_field(field) **CustomElectroMagneticField** -- Arbitrary combined field. The callback must return all six components ``[Bx, By, Bz, Ex, Ey, Ez]``. @@ -130,6 +167,19 @@ Combined magnetic and electric fields. field.field_function = my_EM_field box.add_field(field) +.. warning:: Performance warning + + Same GIL overhead as ``CustomMagneticField``. See :ref:`user_guide_fields_performance`. + +**MappedElectroMagneticField** -- Combined B and E field, each defined on its own independent regular 3D Cartesian grid. The two grids do not need to share the same resolution or spatial extent. + +.. code-block:: python + + field = fields.MappedElectroMagneticField(name="EM_mapped") + field.field_matrix_B = b_matrix # columns: [x, y, z, Bx, By, Bz], values in Geant4 internal units + field.field_matrix_E = e_matrix # columns: [x, y, z, Ex, Ey, Ez], values in Geant4 internal units + box.add_field(field) + Integration accuracy parameters -------------------------------- @@ -244,7 +294,7 @@ Native field types (``UniformMagneticField``, ``UniformElectricField``, ``Quadru Custom fields (``CustomMagneticField``, ``CustomElectricField``, ``CustomElectroMagneticField``) call a Python function for **every evaluation** of ``GetFieldValue`` during tracking. This means that the GIL is acquired and released on every call, which can significantly slow down the simulation, especially in multithreaded mode where all threads serialize through the Python callback. -**Recommendation:** Use native types whenever possible. +**Recommendation:** Use native types whenever possible. For spatially varying fields, use a mapped field type (``MappedMagneticField``, ``MappedElectricField``, ``MappedElectroMagneticField``): define the field on a grid once at setup and let the C++ interpolator handle all evaluations during tracking with zero Python overhead. .. I am keeping this as a comment for now because sources are not serialized yet, so the round-trip for a full simulation object does not work. .. Serialization @@ -271,9 +321,14 @@ The field implementation is covered by the ``test099_fields_*`` tests in ``openg - ``test099_fields_analytical_E`` -- Uniform E field vs analytical energy gain. - ``test099_fields_custom_vs_native_B`` -- Custom trampoline B vs native G4 (bit-identical). - ``test099_fields_custom_vs_native_E`` -- Custom trampoline E vs native G4 (bit-identical). -- ``test099_fields_multi_volume_refresh`` -- Testing field updates when volumes are dynamically modified between runs. -- ``test099_fields_repeated_placements`` -- Testing field behaviour with repeated physical placements of the same logical volume. -- ``test099_fields_rotated_volume`` -- Test that fields rotate with their attached volume. +- ``test099_fields_mapped_vs_uniform_B`` -- MappedMagneticField (constant grid) vs UniformMagneticField. +- ``test099_fields_mapped_vs_uniform_E`` -- MappedElectricField (constant grid) vs UniformElectricField. +- ``test099_fields_multi_volume_refresh`` -- Uniform field shared across two volumes; one is dynamically rotated between runs. +- ``test099_fields_mapped_multi_volume_refresh`` -- Same as above with a MappedMagneticField. +- ``test099_fields_repeated_placements`` -- Uniform field on a single box vs the same total depth split into repeated slabs. +- ``test099_fields_mapped_repeated_placements`` -- Same as above with a MappedMagneticField. +- ``test099_fields_rotated_volume`` -- Uniform field shared between an unrotated and a rotated volume. +- ``test099_fields_mapped_rotated_volume`` -- Same as above with a MappedMagneticField. - ``test099_fields_serialization`` -- Round-trip serialization for all non-custom types. - ``test099_fields_api`` -- API guards. @@ -283,6 +338,7 @@ Class reference .. autoclass:: opengate.geometry.fields.FieldBase :members: + :no-index: .. autoclass:: opengate.geometry.fields.UniformMagneticField :members: @@ -296,11 +352,23 @@ Class reference .. autoclass:: opengate.geometry.fields.CustomMagneticField :members: +.. autoclass:: opengate.geometry.fields.MappedMagneticField + :members: + .. autoclass:: opengate.geometry.fields.UniformElectricField :members: .. autoclass:: opengate.geometry.fields.CustomElectricField :members: +.. autoclass:: opengate.geometry.fields.MappedElectricField + :members: + +.. autoclass:: opengate.geometry.fields.UniformElectroMagneticField + :members: + .. autoclass:: opengate.geometry.fields.CustomElectroMagneticField :members: + +.. autoclass:: opengate.geometry.fields.MappedElectroMagneticField + :members: From 15689937cd5aa696e38ac48c9f7f97abff788152 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Thu, 7 May 2026 08:54:06 +0200 Subject: [PATCH 17/28] Fix test description for MappedElectricField vs UniformElectricField --- .../src/geometry/test099_fields_mapped_vs_uniform_E.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py b/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py index 3931a08e42..36bc488c62 100644 --- a/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py +++ b/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py @@ -1,9 +1,9 @@ #!/usr/bin/env python3 """ -Test 099 - MappedMagneticField with uniform grid vs UniformMagneticField. +Test 099 - MappedElectricField with uniform grid vs UniformElectricField. -Places two boxes side-by-side: one with UniformMagneticField, the other with -a MappedMagneticField whose grid contains a constant B value everywhere. A +Places two boxes side-by-side: one with UniformElectricField, the other with +a MappedElectricField whose grid contains a constant E value everywhere. A proton source fires into each box. Checks that: From fee57e6c019587b0287cc9935220ad214b007b45 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Thu, 7 May 2026 08:58:50 +0200 Subject: [PATCH 18/28] Add round-trip serialization tests for MappedMagneticField, MappedElectricField, and MappedElectroMagneticField --- .../geometry/test099_fields_serialization.py | 71 +++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/opengate/tests/src/geometry/test099_fields_serialization.py b/opengate/tests/src/geometry/test099_fields_serialization.py index 531f3e7f22..92e70b683e 100644 --- a/opengate/tests/src/geometry/test099_fields_serialization.py +++ b/opengate/tests/src/geometry/test099_fields_serialization.py @@ -6,6 +6,10 @@ and that custom fields correctly refuse it. """ +import itertools + +import numpy as np + import opengate as gate from opengate.geometry import fields from opengate.tests import utility @@ -166,4 +170,71 @@ def _make_sim_with_box(): is_ok = is_ok and ok print(f"Multi-volume round-trip: {'OK' if ok else 'FAIL'}") + # Helper: minimal 2×2×2 grid (required for trilinear interpolation) + def _make_grid(half, Fx, Fy, Fz): + corners = list(itertools.product([-half, half], repeat=3)) + return np.array([[x, y, z, Fx, Fy, Fz] for x, y, z in corners]) + + box_half = 25 * g4_cm + + # MappedMagneticField round-trip + sim, box = _make_sim_with_box() + field = fields.MappedMagneticField(name="B_mapped") + field.field_matrix = _make_grid(box_half, 0.0, 1.0 * g4_tesla, 0.0) + field.interpolation = "nearest" + box.add_field(field) + + d = sim.to_dictionary() + sim2 = gate.Simulation() + sim2.from_dictionary(d) + restored = sim2.volume_manager.fields["B_mapped"] + ok = ( + isinstance(restored, fields.MappedMagneticField) + and np.array_equal(restored.field_matrix, field.field_matrix) + and restored.interpolation == "nearest" + and restored.attached_to == ["box"] + ) + is_ok = is_ok and ok + print(f"Mapped B round-trip: {'OK' if ok else 'FAIL'}") + + # MappedElectricField round-trip + sim, box = _make_sim_with_box() + field = fields.MappedElectricField(name="E_mapped") + field.field_matrix = _make_grid(box_half, 0.0, 0.0, 1e6 * g4_volt / g4_m) + box.add_field(field) + + d = sim.to_dictionary() + sim2 = gate.Simulation() + sim2.from_dictionary(d) + restored = sim2.volume_manager.fields["E_mapped"] + ok = ( + isinstance(restored, fields.MappedElectricField) + and np.array_equal(restored.field_matrix, field.field_matrix) + and restored.interpolation == "trilinear" + and restored.attached_to == ["box"] + ) + is_ok = is_ok and ok + print(f"Mapped E round-trip: {'OK' if ok else 'FAIL'}") + + # MappedElectroMagneticField round-trip + sim, box = _make_sim_with_box() + field = fields.MappedElectroMagneticField(name="EB_mapped") + field.field_matrix_B = _make_grid(box_half, 0.0, 2.0 * g4_tesla, 0.0) + field.field_matrix_E = _make_grid(box_half, 0.0, 0.0, 5e5 * g4_volt / g4_m) + box.add_field(field) + + d = sim.to_dictionary() + sim2 = gate.Simulation() + sim2.from_dictionary(d) + restored = sim2.volume_manager.fields["EB_mapped"] + ok = ( + isinstance(restored, fields.MappedElectroMagneticField) + and np.array_equal(restored.field_matrix_B, field.field_matrix_B) + and np.array_equal(restored.field_matrix_E, field.field_matrix_E) + and restored.interpolation == "trilinear" + and restored.attached_to == ["box"] + ) + is_ok = is_ok and ok + print(f"Mapped EB round-trip: {'OK' if ok else 'FAIL'}") + utility.test_ok(is_ok) From 32bc5fe605750b705bcd1e1ee52026a9d89fb544 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Wed, 27 May 2026 11:56:33 +0200 Subject: [PATCH 19/28] Bind new stepper classes --- .../g4_bindings/pyG4BogackiShampine23.cpp | 25 +++++++++++++++++++ .../g4_bindings/pyG4BogackiShampine45.cpp | 25 +++++++++++++++++++ .../g4_bindings/pyG4CashKarpRKF45.cpp | 25 +++++++++++++++++++ .../g4_bindings/pyG4DormandPrince745.cpp | 25 +++++++++++++++++++ .../g4_bindings/pyG4DormandPrinceRK56.cpp | 25 +++++++++++++++++++ .../g4_bindings/pyG4DormandPrinceRK78.cpp | 25 +++++++++++++++++++ .../g4_bindings/pyG4ExactHelixStepper.cpp | 25 +++++++++++++++++++ .../g4_bindings/pyG4NystromRK4.cpp | 25 +++++++++++++++++++ core/opengate_core/opengate_core.cpp | 24 ++++++++++++++++++ 9 files changed, 224 insertions(+) create mode 100644 core/opengate_core/g4_bindings/pyG4BogackiShampine23.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4BogackiShampine45.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4CashKarpRKF45.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4DormandPrince745.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4DormandPrinceRK56.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4DormandPrinceRK78.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4ExactHelixStepper.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4NystromRK4.cpp diff --git a/core/opengate_core/g4_bindings/pyG4BogackiShampine23.cpp b/core/opengate_core/g4_bindings/pyG4BogackiShampine23.cpp new file mode 100644 index 0000000000..6d2e2776d9 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4BogackiShampine23.cpp @@ -0,0 +1,25 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "G4BogackiShampine23.hh" +#include "G4EquationOfMotion.hh" +#include "G4MagIntegratorStepper.hh" + +void init_G4BogackiShampine23(py::module &m) { + // G4BogackiShampine23 inherits from G4MagIntegratorStepper + py::class_>( + m, "G4BogackiShampine23") + + .def(py::init()) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4BogackiShampine45.cpp b/core/opengate_core/g4_bindings/pyG4BogackiShampine45.cpp new file mode 100644 index 0000000000..c18a0bf441 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4BogackiShampine45.cpp @@ -0,0 +1,25 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "G4BogackiShampine45.hh" +#include "G4EquationOfMotion.hh" +#include "G4MagIntegratorStepper.hh" + +void init_G4BogackiShampine45(py::module &m) { + // G4BogackiShampine45 inherits from G4MagIntegratorStepper + py::class_>( + m, "G4BogackiShampine45") + + .def(py::init()) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4CashKarpRKF45.cpp b/core/opengate_core/g4_bindings/pyG4CashKarpRKF45.cpp new file mode 100644 index 0000000000..7525739a30 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4CashKarpRKF45.cpp @@ -0,0 +1,25 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "G4CashKarpRKF45.hh" +#include "G4EquationOfMotion.hh" +#include "G4MagIntegratorStepper.hh" + +void init_G4CashKarpRKF45(py::module &m) { + // G4CashKarpRKF45 inherits from G4MagIntegratorStepper + py::class_>(m, + "G4CashKarpRKF45") + + .def(py::init()) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4DormandPrince745.cpp b/core/opengate_core/g4_bindings/pyG4DormandPrince745.cpp new file mode 100644 index 0000000000..0f4f3611c4 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4DormandPrince745.cpp @@ -0,0 +1,25 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "G4DormandPrince745.hh" +#include "G4EquationOfMotion.hh" +#include "G4MagIntegratorStepper.hh" + +void init_G4DormandPrince745(py::module &m) { + // G4DormandPrince745 inherits from G4MagIntegratorStepper + py::class_>( + m, "G4DormandPrince745") + + .def(py::init()) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4DormandPrinceRK56.cpp b/core/opengate_core/g4_bindings/pyG4DormandPrinceRK56.cpp new file mode 100644 index 0000000000..ae44bb5731 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4DormandPrinceRK56.cpp @@ -0,0 +1,25 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "G4DormandPrinceRK56.hh" +#include "G4EquationOfMotion.hh" +#include "G4MagIntegratorStepper.hh" + +void init_G4DormandPrinceRK56(py::module &m) { + // G4DormandPrinceRK56 inherits from G4MagIntegratorStepper + py::class_>( + m, "G4DormandPrinceRK56") + + .def(py::init()) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4DormandPrinceRK78.cpp b/core/opengate_core/g4_bindings/pyG4DormandPrinceRK78.cpp new file mode 100644 index 0000000000..c3cd9f4b6f --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4DormandPrinceRK78.cpp @@ -0,0 +1,25 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "G4DormandPrinceRK78.hh" +#include "G4EquationOfMotion.hh" +#include "G4MagIntegratorStepper.hh" + +void init_G4DormandPrinceRK78(py::module &m) { + // G4DormandPrinceRK78 inherits from G4MagIntegratorStepper + py::class_>( + m, "G4DormandPrinceRK78") + + .def(py::init()) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4ExactHelixStepper.cpp b/core/opengate_core/g4_bindings/pyG4ExactHelixStepper.cpp new file mode 100644 index 0000000000..4cd83bc253 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4ExactHelixStepper.cpp @@ -0,0 +1,25 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "G4ExactHelixStepper.hh" +#include "G4MagIntegratorStepper.hh" +#include "G4Mag_EqRhs.hh" + +void init_G4ExactHelixStepper(py::module &m) { + // G4ExactHelixStepper <- G4MagHelicalStepper <- G4MagIntegratorStepper. + py::class_>( + m, "G4ExactHelixStepper") + + .def(py::init()) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4NystromRK4.cpp b/core/opengate_core/g4_bindings/pyG4NystromRK4.cpp new file mode 100644 index 0000000000..1add9a97f6 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4NystromRK4.cpp @@ -0,0 +1,25 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "G4MagIntegratorStepper.hh" +#include "G4Mag_EqRhs.hh" +#include "G4NystromRK4.hh" + +void init_G4NystromRK4(py::module &m) { + // G4NystromRK4 inherits from G4MagIntegratorStepper + py::class_>(m, "G4NystromRK4") + + .def(py::init(), py::arg("equation"), + py::arg("distanceConstField") = 0.0) + + ; +} diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index e9a6c1ec1a..f11b571323 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -206,6 +206,22 @@ void init_G4MagErrorStepper(py::module &); void init_G4ClassicalRK4(py::module &); +void init_G4DormandPrince745(py::module &); + +void init_G4DormandPrinceRK56(py::module &); + +void init_G4DormandPrinceRK78(py::module &); + +void init_G4BogackiShampine23(py::module &); + +void init_G4BogackiShampine45(py::module &); + +void init_G4CashKarpRKF45(py::module &); + +void init_G4NystromRK4(py::module &); + +void init_G4ExactHelixStepper(py::module &); + void init_G4VIntegrationDriver(py::module &); void init_G4MagInt_Driver(py::module &); @@ -592,6 +608,14 @@ PYBIND11_MODULE(opengate_core, m) { init_G4MagIntegratorStepper(m); init_G4MagErrorStepper(m); init_G4ClassicalRK4(m); + init_G4DormandPrince745(m); + init_G4DormandPrinceRK56(m); + init_G4DormandPrinceRK78(m); + init_G4BogackiShampine23(m); + init_G4BogackiShampine45(m); + init_G4CashKarpRKF45(m); + init_G4NystromRK4(m); + init_G4ExactHelixStepper(m); init_G4VIntegrationDriver(m); init_G4MagInt_Driver(m); init_G4ChordFinder(m); From a80a74c73cfad0c71ed0345cd4d0173bf4d71f61 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Wed, 27 May 2026 12:04:46 +0200 Subject: [PATCH 20/28] Add stepper type parameter to FieldBase class + guards for magnetic-only steppers --- opengate/geometry/fields.py | 84 ++++++++++++++++++++++++++++--------- 1 file changed, 65 insertions(+), 19 deletions(-) diff --git a/opengate/geometry/fields.py b/opengate/geometry/fields.py index 04cd6821b9..1bba79cf65 100644 --- a/opengate/geometry/fields.py +++ b/opengate/geometry/fields.py @@ -28,6 +28,7 @@ class FieldBase(GateObject): """Base class for electric and magnetic fields.""" # hints for IDE + stepper: str step_minimum: float delta_chord: float delta_one_step: float @@ -36,6 +37,18 @@ class FieldBase(GateObject): max_epsilon_step: float user_info_defaults = { + "stepper": ( + "DormandPrince745", + { + "doc": ( + "Integration stepper type. " + "General-purpose (any field): 'DormandPrince745' (default), 'CashKarpRKF45', " + "'BogackiShampine45', 'BogackiShampine23', 'DormandPrinceRK56', " + "'DormandPrinceRK78', 'ClassicalRK4'. " + "Magnetic-only: 'NystromRK4', 'ExactHelixStepper'." + ), + }, + ), "step_minimum": ( 1e-2 * g4_units.mm, { @@ -143,29 +156,28 @@ def _make_g4_transforms(self): [rot_np_as_g4(r) for r in rotations_np], ) - @staticmethod - def _validate_field_function(func, class_name, n_components): - """Validate that func is callable and returns exactly n_components values.""" - if func is None: - raise ValueError(f"field_function must be provided for {class_name}") - if not callable(func): - raise TypeError("field_function must be a callable function") - result = list(func(0, 0, 0, 0)) - if len(result) != n_components: + def _validate_stepper(self) -> None: + """Raise ValueError if the current stepper is incompatible with this field type.""" + if self.stepper not in _stepper_map: + raise ValueError( + f"Unknown stepper '{self.stepper}'. " + f"Choose from {list(_stepper_map.keys())}." + ) + if self.stepper in _magnetic_only_steppers and self.field_changes_energy: raise ValueError( - f"{class_name}: field_function must return {n_components} components, " - f"got {len(result)}" + f"Stepper '{self.stepper}' only supports pure magnetic fields. " + "Choose a general-purpose stepper for electric or EM fields." ) def _build_field_manager( self, inner_field, gate_field, equation_cls, n_vars, volume_obj ): """Build equation/stepper/driver/chord_finder/fm, record runtime objects, return fm.""" + self._validate_stepper() self.g4_field = gate_field self.g4_equation_of_motion = equation_cls(gate_field) - self.g4_integrator_stepper = g4.G4ClassicalRK4( - self.g4_equation_of_motion, n_vars - ) + stepper_factory = _stepper_map[self.stepper] + self.g4_integrator_stepper = stepper_factory(self.g4_equation_of_motion, n_vars) self.g4_integration_driver = g4.G4MagInt_Driver( self.step_minimum, self.g4_integrator_stepper, n_vars, 0 ) @@ -328,7 +340,7 @@ def _create_inner_field(self): the attached volume and must return [Bx, By, Bz] in that same local frame. The base class rotates the result to world coordinates. """ - self._validate_field_function(self.field_function, "CustomMagneticField", 3) + _validate_field_function(self.field_function, "CustomMagneticField", 3) class _PyMagneticField(g4.G4MagneticField): def __init__(inner_self, callback): @@ -422,7 +434,7 @@ class CustomElectricField(ElectricField): } def _create_inner_field(self): - self._validate_field_function(self.field_function, "CustomElectricField", 3) + _validate_field_function(self.field_function, "CustomElectricField", 3) class _PyElectricField(g4.G4ElectricField): def __init__(inner_self, callback): @@ -488,9 +500,7 @@ class CustomElectroMagneticField(ElectroMagneticField): } def _create_inner_field(self): - self._validate_field_function( - self.field_function, "CustomElectroMagneticField", 6 - ) + _validate_field_function(self.field_function, "CustomElectroMagneticField", 6) class _PyElectroMagneticField(g4.G4ElectroMagneticField): def __init__(inner_self, callback): @@ -511,6 +521,7 @@ def to_dictionary(self): ) +# Helper function to parse field_matrix for mapped fields and extract grid parameters and field value arrays def _parse_field_matrix(mat, class_name): """Parse a field matrix into grid parameters and sorted field value arrays. @@ -566,6 +577,41 @@ def _parse_field_matrix(mat, class_name): return nx, ny, nz, x0, y0, z0, dx, dy, dz, field_cols +# Helper function to validate user-provided field_function for custom fields +def _validate_field_function( + func: Any, + class_name: str, + n_components: int, +) -> None: + """Validate that field_function is a callable that returns the expected number of components.""" + if func is None: + raise ValueError(f"field_function must be provided for {class_name}") + if not callable(func): + raise TypeError("field_function must be a callable function") + result = list(func(0, 0, 0, 0)) + if len(result) != n_components: + raise ValueError( + f"{class_name}: field_function must return {n_components} components, " + f"got {len(result)}" + ) + + +# Steppers that only work with pure magnetic fields +_magnetic_only_steppers = {"NystromRK4", "ExactHelixStepper"} + +# Stepper factory +_stepper_map = { + "DormandPrince745": lambda eq, n: g4.G4DormandPrince745(eq, n), + "CashKarpRKF45": lambda eq, n: g4.G4CashKarpRKF45(eq, n), + "BogackiShampine45": lambda eq, n: g4.G4BogackiShampine45(eq, n), + "BogackiShampine23": lambda eq, n: g4.G4BogackiShampine23(eq, n), + "DormandPrinceRK56": lambda eq, n: g4.G4DormandPrinceRK56(eq, n), + "DormandPrinceRK78": lambda eq, n: g4.G4DormandPrinceRK78(eq, n), + "ClassicalRK4": lambda eq, n: g4.G4ClassicalRK4(eq, n), + "NystromRK4": lambda eq, n: g4.G4NystromRK4(eq), + "ExactHelixStepper": lambda eq, n: g4.G4ExactHelixStepper(eq), +} + _interp_map = { "trilinear": g4.GateGridInterpolationMethod.Trilinear, "nearest": g4.GateGridInterpolationMethod.Nearest, From 780d92e4c01999053b5cdb8269e66d7df9f8f777 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Wed, 27 May 2026 13:45:23 +0200 Subject: [PATCH 21/28] Update TODOs in fields.py --- opengate/geometry/fields.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/opengate/geometry/fields.py b/opengate/geometry/fields.py index 1bba79cf65..d490fca4a2 100644 --- a/opengate/geometry/fields.py +++ b/opengate/geometry/fields.py @@ -14,10 +14,6 @@ from ..utility import g4_units # ! ======= KNOWN TODO'S ======== -# ! - implement the possibility of choosing the stepper type and equation type -# ! WIP - implement mapped fields (e.g., from a CSV file) -# ! - Overhead for custom fields implementation: every GetFieldValue call crosses c++ -> Python -> c++, acquiring the GIL each time, -# ! which is very inefficient. Need to implement a more efficient way on the C++ side. # ! - in MT mode, create_field_manager() is called per thread and overwrites shared instance attrs (g4_field, etc.), # ! so any code that relies on those attributes being available later on will get an arbitrary thread's copy. # ! Not sure if this is really an issue or not, but worth investigating. @@ -42,9 +38,9 @@ class FieldBase(GateObject): { "doc": ( "Integration stepper type. " - "General-purpose (any field): 'DormandPrince745' (default), 'CashKarpRKF45', " + "General-purpose (any field): 'DormandPrince745' (default), 'ClassicalRK4', 'CashKarpRKF45', " "'BogackiShampine45', 'BogackiShampine23', 'DormandPrinceRK56', " - "'DormandPrinceRK78', 'ClassicalRK4'. " + "'DormandPrinceRK78'. " "Magnetic-only: 'NystromRK4', 'ExactHelixStepper'." ), }, From b19b7844736226ae6e12af31fe26349b27b0e88d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Wed, 27 May 2026 13:50:42 +0200 Subject: [PATCH 22/28] Update documentation --- docs/source/user_guide/user_guide_fields.rst | 32 +++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/docs/source/user_guide/user_guide_fields.rst b/docs/source/user_guide/user_guide_fields.rst index 5081909b55..c3fd1ed7c4 100644 --- a/docs/source/user_guide/user_guide_fields.rst +++ b/docs/source/user_guide/user_guide_fields.rst @@ -181,7 +181,7 @@ Combined magnetic and electric fields. box.add_field(field) -Integration accuracy parameters +Integration and accuracy parameters -------------------------------- All field types inherit the following parameters that control the numerical integration of the equation of motion. The defaults are suitable for most cases, but they can be tuned for better accuracy or performance. @@ -193,6 +193,9 @@ All field types inherit the following parameters that control the numerical inte * - Parameter - Default - Description + * - ``stepper`` + - ``DormandPrince745`` + - Stepping algorithm used for integrating the equation of motion. See :ref:`steppers` for available options. * - ``step_minimum`` - 0.01 mm - Minimum step size for the chord finder. @@ -218,11 +221,36 @@ Example: field = fields.UniformMagneticField(name="B") field.field_vector = [0, 0, 1 * tesla] + field.stepper = "DormandPrince745" field.delta_chord = 0.01 * mm field.step_minimum = 0.1 * mm box.add_field(field) +.. _steppers: + +Steppers +~~~~~~~~ + +For general fields (electric and/or magnetic components): + +- ``DormandPrince745`` (default) +- ``ClassicalRK4`` +- ``CashKarpRKF45`` +- ``BogackiShampine45`` +- ``BogackiShampine23`` +- ``DormandPrinceRK56`` +- ``DormandPrinceRK78`` + +For purely magnetic fields, the following additional steppers are available: + +- ``NystromRK4`` +- ``ExactHelicalStepper`` + +Please refer to the `Geant4 documentation `__ for details on the characteristics and recommended use cases for each stepper type. + + + Attaching fields to volumes ---------------------------- @@ -331,6 +359,8 @@ The field implementation is covered by the ``test099_fields_*`` tests in ``openg - ``test099_fields_mapped_rotated_volume`` -- Same as above with a MappedMagneticField. - ``test099_fields_serialization`` -- Round-trip serialization for all non-custom types. - ``test099_fields_api`` -- API guards. +- ``test099_fields_stepper_em`` -- Comparison of different steppers for a uniform E field. +- ``test099_fields_stepper`` -- Comparison of different steppers for a uniform B field. Class reference From f418ae266480936eb0994b4e2cb6323eedf6c4a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Wed, 27 May 2026 15:16:31 +0200 Subject: [PATCH 23/28] Add tests --- .../src/geometry/test099_fields_stepper.py | 149 +++++++++++++++++ .../src/geometry/test099_fields_stepper_em.py | 157 ++++++++++++++++++ 2 files changed, 306 insertions(+) create mode 100644 opengate/tests/src/geometry/test099_fields_stepper.py create mode 100644 opengate/tests/src/geometry/test099_fields_stepper_em.py diff --git a/opengate/tests/src/geometry/test099_fields_stepper.py b/opengate/tests/src/geometry/test099_fields_stepper.py new file mode 100644 index 0000000000..ed3ff55766 --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_stepper.py @@ -0,0 +1,149 @@ +#!/usr/bin/env python3 +""" +Test 099 - Stepper selection: all implemented steppers with a uniform B field. + +Places one box per stepper side-by-side in a single simulation, each with the +same uniform magnetic field. A proton source fires into each box. + +Covers all general-purpose and magnetic-only steppers (NystromRK4, +ExactHelixStepper). Checks that: + - Every stepper agrees with the analytical cyclotron radius. + - Every stepper agrees with DormandPrince745 (the reference) in exit + position and kinetic energy. + +See test099_fields_stepper_em.py for non-pure-B field coverage and +magnetic-only stepper guard tests. +""" + +import numpy as np +import uproot + +import opengate as gate +from opengate.geometry import fields +from opengate.tests import utility + +from test099_fields_helpers import ( + g4_m, + g4_cm, + g4_mm, + g4_tesla, + g4_MeV, + g4_eplus, + PROTON_MASS, + cyclotron_radius, +) + +REF_STEPPER = "DormandPrince745" +ALL_STEPPERS = list(fields._stepper_map.keys()) + +VISU = False + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, output_folder="test099_fields") + + By = 5 * g4_tesla + T = 10 * g4_MeV + + n = len(ALL_STEPPERS) + spacing = 80 * g4_cm + + x_positions = [(i - n // 2) * spacing for i in range(n)] + + sim = gate.Simulation() + sim.g4_verbose = False + sim.visu = False + sim.number_of_threads = 1 + sim.random_seed = 42 + sim.output_dir = paths.output + + if VISU: + sim.visu = True + sim.visu_type = "qt" + sim.visu_commands.append("/vis/scene/endOfEventAction accumulate") + sim.visu_commands.append("/vis/scene/add/trajectories smooth") + sim.visu_commands.append("/vis/scene/add/magneticField 20 fullArrow") + + world = sim.world + world.size = [(n + 1) * spacing, 1 * g4_m, 1 * g4_m] + world.material = "G4_Galactic" + + for stepper_name, x in zip(ALL_STEPPERS, x_positions): + box = sim.add_volume("BoxVolume", f"box_{stepper_name}") + box.size = [50 * g4_cm, 50 * g4_cm, 50 * g4_cm] + box.translation = [x, 0, 0] + box.material = "G4_Galactic" + + field = fields.UniformMagneticField(name=f"B_{stepper_name}") + field.field_vector = [0, By, 0] + field.stepper = stepper_name + box.add_field(field) + + src = sim.add_source("GenericSource", f"src_{stepper_name}") + src.particle = "proton" + src.n = 1 + src.energy.type = "mono" + src.energy.mono = T + src.position.type = "point" + src.position.translation = [x, 0, -100 * g4_cm] + src.direction.type = "momentum" + src.direction.momentum = [0, 0, 1] + + phsp = sim.add_actor("PhaseSpaceActor", f"phsp_{stepper_name}") + phsp.attached_to = f"box_{stepper_name}" + phsp.attributes = ["PostKineticEnergy", "PostPosition"] + phsp.output_filename = paths.output / f"test099_stepper_{stepper_name}.root" + phsp.steps_to_store = "exiting" + + sim.run() + + # --- Collect results --- + results = {} + for stepper_name, x in zip(ALL_STEPPERS, x_positions): + df = uproot.open(str(paths.output / f"test099_stepper_{stepper_name}.root"))[ + f"phsp_{stepper_name};1" + ].arrays(library="pd") + results[stepper_name] = ( + df["PostPosition_X"].values - x, # shift to box-local frame + df["PostPosition_Z"].values, + df["PostKineticEnergy"].values, + ) + + # --- Analytical reference --- + r = cyclotron_radius(T, By, PROTON_MASS, 1 * g4_eplus) + box_half_z = 250 * g4_mm + cx, cz = 0.0 - r, -box_half_z + + r_TOL = 1e-3 * g4_mm + e_TOL = 1e-3 * g4_MeV + + is_ok = True + + # Each stepper vs. analytical + for stepper_name, (pos_x, pos_z, KE) in results.items(): + residual = np.sqrt((pos_x - cx) ** 2 + (pos_z - cz) ** 2) - r + ok_pos = bool(np.all(np.abs(residual) < r_TOL)) + ok_e = bool(np.all(np.abs(KE - T) < e_TOL)) + print( + f"{stepper_name:25s} vs analytical: " + f"residual={residual[0] / g4_mm:.6f} mm " + f"OK_pos={ok_pos} OK_E={ok_e}" + ) + is_ok = is_ok and ok_pos and ok_e + + # Each stepper vs. reference (DormandPrince745) + ref_x, ref_z, ref_KE = results[REF_STEPPER] + for stepper_name, (pos_x, pos_z, KE) in results.items(): + if stepper_name == REF_STEPPER: + continue + dx = np.abs(pos_x - ref_x).max() + dz = np.abs(pos_z - ref_z).max() + dKE = np.abs(KE - ref_KE).max() + ok_cross = bool(dx < r_TOL and dz < r_TOL and dKE < e_TOL) + print( + f"{stepper_name:25s} vs {REF_STEPPER}: " + f"dx={dx / g4_mm:.6f} mm dz={dz / g4_mm:.6f} mm " + f"dKE={dKE:.6f} MeV OK={ok_cross}" + ) + is_ok = is_ok and ok_cross + + utility.test_ok(is_ok) diff --git a/opengate/tests/src/geometry/test099_fields_stepper_em.py b/opengate/tests/src/geometry/test099_fields_stepper_em.py new file mode 100644 index 0000000000..1484cb4ced --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_stepper_em.py @@ -0,0 +1,157 @@ +#!/usr/bin/env python3 +""" +Test 099 - Stepper selection with non-pure-B fields. + +Test 1 All general-purpose steppers with a uniform transverse electric field + (Ex along X). All steppers must agree with each other in exit position + and kinetic energy (cross-stepper consistency). Physical correctness of + the field implementation is validated separately in test099_fields_analytical_E. + +Test 2 Magnetic-only steppers (NystromRK4, ExactHelixStepper) paired with + UniformElectricField or UniformElectroMagneticField must raise a + ValueError. +""" + +import numpy as np +import uproot + +import opengate as gate +from opengate.geometry import fields +from opengate.tests import utility + +from test099_fields_helpers import ( + g4_m, + g4_cm, + g4_mm, + g4_MeV, + g4_volt, + g4_tesla, + PROTON_MASS, +) + +GENERAL_PURPOSE = [ + s for s in fields._stepper_map if s not in fields._magnetic_only_steppers +] +MAGNETIC_ONLY = sorted(fields._magnetic_only_steppers) +REF_STEPPER = "DormandPrince745" + +VISU = True + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, output_folder="test099_fields") + + T = 100 * g4_MeV + L = 500 * g4_mm + Ex = 1e6 * g4_volt / g4_cm + By = 3 * g4_tesla + + # Test 1: all general-purpose steppers with transverse E field + n = len(GENERAL_PURPOSE) + spacing = 80 * g4_cm + x_positions = [(i - n // 2) * spacing for i in range(n)] + + sim = gate.Simulation() + sim.g4_verbose = False + sim.visu = False + sim.number_of_threads = 1 + sim.random_seed = 42 + sim.output_dir = paths.output + + if VISU: + sim.visu = True + sim.visu_type = "qt" + sim.visu_commands.append("/vis/scene/endOfEventAction accumulate") + sim.visu_commands.append("/vis/scene/add/trajectories smooth") + sim.visu_commands.append("/vis/scene/add/electricField 20 fullArrow") + + sim.world.size = [(n + 1) * spacing, 1 * g4_m, 1 * g4_m] + sim.world.material = "G4_Galactic" + + for stepper_name, x in zip(GENERAL_PURPOSE, x_positions): + box = sim.add_volume("BoxVolume", f"box_{stepper_name}") + box.size = [50 * g4_cm, 50 * g4_cm, 50 * g4_cm] + box.translation = [x, 0, 0] + box.material = "G4_Galactic" + + field = fields.UniformElectricField(name=f"E_{stepper_name}") + field.field_vector = [Ex, 0, 0] + field.stepper = stepper_name + box.add_field(field) + + src = sim.add_source("GenericSource", f"src_{stepper_name}") + src.particle = "proton" + src.n = 1 + src.energy.type = "mono" + src.energy.mono = T + src.position.type = "point" + src.position.translation = [x, 0, -100 * g4_cm] + src.direction.type = "momentum" + src.direction.momentum = [0, 0, 1] + + phsp = sim.add_actor("PhaseSpaceActor", f"phsp_{stepper_name}") + phsp.attached_to = f"box_{stepper_name}" + phsp.attributes = ["PostKineticEnergy", "PostPosition"] + phsp.output_filename = paths.output / f"test099_stepper_em_{stepper_name}.root" + phsp.steps_to_store = "exiting" + + sim.run() + + x_TOL = 0.01 * g4_mm + e_TOL = 1e-3 * g4_MeV + is_ok = True + results = {} + + for stepper_name, x in zip(GENERAL_PURPOSE, x_positions): + df = uproot.open(str(paths.output / f"test099_stepper_em_{stepper_name}.root"))[ + f"phsp_{stepper_name};1" + ].arrays(library="pd") + x_exit = float(df["PostPosition_X"].values[0]) - x + KE_exit = float(df["PostKineticEnergy"].values[0]) + results[stepper_name] = (x_exit, KE_exit) + print( + f"{stepper_name:25s} x={x_exit / g4_mm:.4f} mm KE={KE_exit / g4_MeV:.4f} MeV" + ) + + ref_x, ref_KE = results[REF_STEPPER] + for stepper_name, (x_exit, KE_exit) in results.items(): + if stepper_name == REF_STEPPER: + continue + dx = abs(x_exit - ref_x) + dKE = abs(KE_exit - ref_KE) + ok_cross = dx < x_TOL and dKE < e_TOL + print( + f"{stepper_name:25s} vs {REF_STEPPER}: " + f"dx={dx / g4_mm:.6f} mm dKE={dKE / g4_MeV:.6f} MeV OK={ok_cross}" + ) + is_ok = is_ok and ok_cross + + # Test 2: magnetic-only steppers must raise ValueError with E/EM fields + guard_cases = [ + ("E only", fields.UniformElectricField, {"field_vector": [0, 0, Ex]}), + ( + "E + B", + fields.UniformElectroMagneticField, + {"field_vector_E": [0, 0, Ex], "field_vector_B": [0, By, 0]}, + ), + ] + + for stepper_name in MAGNETIC_ONLY: + for label, FieldClass, field_kwargs in guard_cases: + fld = FieldClass(name="guard_field") + for attr, val in field_kwargs.items(): + setattr(fld, attr, val) + fld.stepper = stepper_name + + try: + fld._validate_stepper() + ok_guard = False + except ValueError as exc: + ok_guard = "only supports pure magnetic fields" in str(exc) + + print( + f"{stepper_name:25s} + {label:6s}: " + f"ValueError raised correctly: {ok_guard}" + ) + is_ok = is_ok and ok_guard + + utility.test_ok(is_ok) From 171727ee00715f17b45915a9d16f5c9d5c290fc2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Wed, 27 May 2026 15:16:31 +0200 Subject: [PATCH 24/28] Add tests --- opengate/tests/src/geometry/test099_fields_stepper_em.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opengate/tests/src/geometry/test099_fields_stepper_em.py b/opengate/tests/src/geometry/test099_fields_stepper_em.py index 1484cb4ced..1429bf94e5 100644 --- a/opengate/tests/src/geometry/test099_fields_stepper_em.py +++ b/opengate/tests/src/geometry/test099_fields_stepper_em.py @@ -35,7 +35,7 @@ MAGNETIC_ONLY = sorted(fields._magnetic_only_steppers) REF_STEPPER = "DormandPrince745" -VISU = True +VISU = False if __name__ == "__main__": paths = utility.get_default_test_paths(__file__, output_folder="test099_fields") From 50fdead90a6bb140f11819875bccaf30921c912b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Sun, 31 May 2026 08:58:53 +0200 Subject: [PATCH 25/28] Trigger CI tests From 909c04fb5f3df805d5f72dfe3bc94b3376bc3ec3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Mon, 1 Jun 2026 07:56:57 +0200 Subject: [PATCH 26/28] Re-trigger test suite From bad550a1b6db2594f859dd2a93168ccd0479d61f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Tue, 2 Jun 2026 08:06:30 +0200 Subject: [PATCH 27/28] Re-trigger test suite From 919b68bbfc0f77693021e8c68d2414545f704140 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc=20Ballestero=20Rib=C3=B3?= Date: Thu, 4 Jun 2026 10:04:21 +0200 Subject: [PATCH 28/28] Relax tolerances in test --- opengate/tests/src/geometry/test099_fields_stepper.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opengate/tests/src/geometry/test099_fields_stepper.py b/opengate/tests/src/geometry/test099_fields_stepper.py index ed3ff55766..95a57749da 100644 --- a/opengate/tests/src/geometry/test099_fields_stepper.py +++ b/opengate/tests/src/geometry/test099_fields_stepper.py @@ -113,8 +113,8 @@ box_half_z = 250 * g4_mm cx, cz = 0.0 - r, -box_half_z - r_TOL = 1e-3 * g4_mm - e_TOL = 1e-3 * g4_MeV + r_TOL = 1e-2 * g4_mm + e_TOL = 1e-2 * g4_MeV is_ok = True