From 55753de4a7a115cef2503a7879675af4c879bc81 Mon Sep 17 00:00:00 2001 From: srmarcballestero Date: Tue, 21 Apr 2026 16:32:41 +0200 Subject: [PATCH] Implement electric and magnetic fields in GATE Co-authored-by: srmarcballestero Co-authored-by: clemencegentner --- .../g4_bindings/pyG4ChordFinder.cpp | 26 + .../g4_bindings/pyG4ClassicalRK4.cpp | 24 + .../g4_bindings/pyG4ElectricField.cpp | 75 +++ .../g4_bindings/pyG4ElectroMagneticField.cpp | 80 +++ .../g4_bindings/pyG4EqMagElectricField.cpp | 28 + .../g4_bindings/pyG4EquationOfMotion.cpp | 23 + core/opengate_core/g4_bindings/pyG4Field.cpp | 80 +++ .../g4_bindings/pyG4FieldManager.cpp | 66 +++ .../g4_bindings/pyG4MagErrorStepper.cpp | 22 + .../g4_bindings/pyG4MagInt_Driver.cpp | 27 + .../g4_bindings/pyG4MagIntegratorStepper.cpp | 28 + .../g4_bindings/pyG4Mag_EqRhs.cpp | 21 + .../g4_bindings/pyG4Mag_UsualEqRhs.cpp | 29 + .../g4_bindings/pyG4MagneticField.cpp | 75 +++ .../g4_bindings/pyG4QuadrupoleMagField.cpp | 24 + .../g4_bindings/pyG4UniformElectricField.cpp | 25 + .../g4_bindings/pyG4UniformMagField.cpp | 28 + .../g4_bindings/pyG4VIntegrationDriver.cpp | 27 + core/opengate_core/opengate_core.cpp | 55 ++ docs/source/user_guide/index.rst | 2 + docs/source/user_guide/user_guide_fields.rst | 272 +++++++++ .../user_guide_reference_fields.rst | 25 + opengate/__init__.py | 1 + opengate/engines.py | 8 + opengate/geometry/fields.py | 553 ++++++++++++++++++ opengate/geometry/volumes.py | 31 + opengate/managers.py | 16 + .../geometry/test099_fields_analytical_B.py | 80 +++ .../geometry/test099_fields_analytical_E.py | 71 +++ .../tests/src/geometry/test099_fields_api.py | 149 +++++ .../test099_fields_custom_vs_native_B.py | 158 +++++ .../test099_fields_custom_vs_native_E.py | 156 +++++ .../src/geometry/test099_fields_helpers.py | 84 +++ .../geometry/test099_fields_serialization.py | 151 +++++ 34 files changed, 2520 insertions(+) create mode 100644 core/opengate_core/g4_bindings/pyG4ChordFinder.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4ClassicalRK4.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4ElectricField.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4ElectroMagneticField.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4EqMagElectricField.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4EquationOfMotion.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4Field.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4FieldManager.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4MagErrorStepper.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4MagInt_Driver.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4MagIntegratorStepper.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4Mag_EqRhs.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4Mag_UsualEqRhs.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4MagneticField.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4QuadrupoleMagField.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4UniformElectricField.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4UniformMagField.cpp create mode 100644 core/opengate_core/g4_bindings/pyG4VIntegrationDriver.cpp create mode 100644 docs/source/user_guide/user_guide_fields.rst create mode 100644 docs/source/user_guide/user_guide_reference_fields.rst create mode 100644 opengate/geometry/fields.py create mode 100644 opengate/tests/src/geometry/test099_fields_analytical_B.py create mode 100644 opengate/tests/src/geometry/test099_fields_analytical_E.py create mode 100644 opengate/tests/src/geometry/test099_fields_api.py create mode 100644 opengate/tests/src/geometry/test099_fields_custom_vs_native_B.py create mode 100644 opengate/tests/src/geometry/test099_fields_custom_vs_native_E.py create mode 100644 opengate/tests/src/geometry/test099_fields_helpers.py create mode 100644 opengate/tests/src/geometry/test099_fields_serialization.py diff --git a/core/opengate_core/g4_bindings/pyG4ChordFinder.cpp b/core/opengate_core/g4_bindings/pyG4ChordFinder.cpp new file mode 100644 index 000000000..d0ea8f7f9 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4ChordFinder.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 + +namespace py = pybind11; + +#include "G4ChordFinder.hh" +#include "G4MagIntegratorStepper.hh" +#include "G4MagneticField.hh" +#include "G4VIntegrationDriver.hh" + +void init_G4ChordFinder(py::module &m) { + py::class_>( + m, "G4ChordFinder") + + .def(py::init()) + .def(py::init()) + + .def("SetDeltaChord", &G4ChordFinder::SetDeltaChord); +} diff --git a/core/opengate_core/g4_bindings/pyG4ClassicalRK4.cpp b/core/opengate_core/g4_bindings/pyG4ClassicalRK4.cpp new file mode 100644 index 000000000..8be387b3c --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4ClassicalRK4.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 "G4ClassicalRK4.hh" +#include "G4EquationOfMotion.hh" +#include "G4MagErrorStepper.hh" + +void init_G4ClassicalRK4(py::module &m) { + // G4ClassicalRK4 inherits from G4MagErrorStepper + py::class_>(m, "G4ClassicalRK4") + + .def(py::init()) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4ElectricField.cpp b/core/opengate_core/g4_bindings/pyG4ElectricField.cpp new file mode 100644 index 000000000..85b91fca2 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4ElectricField.cpp @@ -0,0 +1,75 @@ +/* -------------------------------------------------- + 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 "G4ElectricField.hh" +#include "G4ElectroMagneticField.hh" + +// Trampoline class to allow Python to override GetFieldValue for electric +// fields. Inherits from G4ElectricField which already implements +// DoesFieldChangeEnergy() = true. +class PyG4ElectricField : public G4ElectricField { +public: + using G4ElectricField::G4ElectricField; + + void GetFieldValue(const G4double Point[4], G4double *field) const override { + // Always initialize output to a safe value: [Bx, By, Bz, Ex, Ey, Ez]. + field[0] = 0.; + field[1] = 0.; + field[2] = 0.; + field[3] = 0.; + field[4] = 0.; + field[5] = 0.; + + py::gil_scoped_acquire gil; + + // Convert Point to Python list + py::list pyPoint; + pyPoint.append(Point[0]); + pyPoint.append(Point[1]); + pyPoint.append(Point[2]); + pyPoint.append(Point[3]); + + // Get the Python override + py::function override = py::get_override(this, "GetFieldValue"); + if (override) { + py::object result = override(pyPoint); + + if (!result.is_none()) { + py::sequence field_seq = result.cast(); + + if (py::len(field_seq) != 3) { + throw std::invalid_argument( + "GetFieldValue for G4ElectricField must return exactly 3 " + "components [Ex, Ey, Ez]"); + } + + // User returned [Ex, Ey, Ez] + field[3] = field_seq[0].cast(); + field[4] = field_seq[1].cast(); + field[5] = field_seq[2].cast(); + } + } + } +}; + +void init_G4ElectricField(py::module &m) { + + py::class_>(m, + "G4ElectricField") + + .def(py::init<>()) + + .def("DoesFieldChangeEnergy", &G4ElectricField::DoesFieldChangeEnergy) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4ElectroMagneticField.cpp b/core/opengate_core/g4_bindings/pyG4ElectroMagneticField.cpp new file mode 100644 index 000000000..fb30746b2 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4ElectroMagneticField.cpp @@ -0,0 +1,80 @@ +/* -------------------------------------------------- + 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 "G4Field.hh" + +// Trampoline class to allow Python to override GetFieldValue for +// electromagnetic fields. +class PyG4ElectroMagneticField : public G4ElectroMagneticField { +public: + using G4ElectroMagneticField::G4ElectroMagneticField; + + void GetFieldValue(const G4double Point[4], G4double *field) const override { + // Always initialize output to a safe value: [Bx, By, Bz, Ex, Ey, Ez]. + field[0] = 0.; + field[1] = 0.; + field[2] = 0.; + field[3] = 0.; + field[4] = 0.; + field[5] = 0.; + + py::gil_scoped_acquire gil; + + // Convert Point to Python list + py::list pyPoint; + pyPoint.append(Point[0]); + pyPoint.append(Point[1]); + pyPoint.append(Point[2]); + pyPoint.append(Point[3]); + + // Try to get the Python override + py::function override = py::get_override(this, "GetFieldValue"); + if (override) { + // Call Python implementation and expect [Bx, By, Bz, Ex, Ey, Ez] + py::object result = override(pyPoint); + if (!result.is_none()) { + py::sequence field_seq = result.cast(); + size_t n = py::len(field_seq); + + if (n != 6) { + throw std::invalid_argument( + "GetFieldValue for G4ElectroMagneticField must return exactly 6 " + "components [Bx, By, Bz, Ex, Ey, Ez]"); + } + + for (size_t i = 0; i < n && i < 6; ++i) { + field[i] = field_seq[i].cast(); + } + } + } + } + + G4bool DoesFieldChangeEnergy() const override { + PYBIND11_OVERRIDE_PURE(G4bool, G4ElectroMagneticField, + DoesFieldChangeEnergy); + } +}; + +void init_G4ElectroMagneticField(py::module &m) { + + py::class_>( + m, "G4ElectroMagneticField") + + .def(py::init<>()) + + .def("DoesFieldChangeEnergy", + &G4ElectroMagneticField::DoesFieldChangeEnergy) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4EqMagElectricField.cpp b/core/opengate_core/g4_bindings/pyG4EqMagElectricField.cpp new file mode 100644 index 000000000..ef849e565 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4EqMagElectricField.cpp @@ -0,0 +1,28 @@ +/* -------------------------------------------------- + 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 "G4ElectroMagneticField.hh" +#include "G4EqMagElectricField.hh" +#include "G4EquationOfMotion.hh" + +void init_G4EqMagElectricField(py::module &m) { + + py::class_>( + m, "G4EqMagElectricField") + + .def(py::init()) + + .def("SetChargeMomentumMass", + &G4EqMagElectricField::SetChargeMomentumMass) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4EquationOfMotion.cpp b/core/opengate_core/g4_bindings/pyG4EquationOfMotion.cpp new file mode 100644 index 000000000..2459bd64e --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4EquationOfMotion.cpp @@ -0,0 +1,23 @@ +/* -------------------------------------------------- + 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 "G4EquationOfMotion.hh" + +void init_G4EquationOfMotion(py::module &m) { + + py::class_>( + m, "G4EquationOfMotion") + + .def("GetFieldObj", + py::overload_cast<>(&G4EquationOfMotion::GetFieldObj, py::const_), + py::return_value_policy::reference_internal); +} diff --git a/core/opengate_core/g4_bindings/pyG4Field.cpp b/core/opengate_core/g4_bindings/pyG4Field.cpp new file mode 100644 index 000000000..e8b8b44a0 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4Field.cpp @@ -0,0 +1,80 @@ +/* -------------------------------------------------- + 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 "G4Field.hh" + +/* + * Trampoline class to allow Python to override GetFieldValue and + * DoesFieldChangeEnergy. This enables users to define custom fields in Python. + * + * The GetFieldValue method receives a point (x, y, z, t) and must return + * the field components at that point. The number of components depends on + * the field type: + * - Magnetic field: 3 components (Bx, By, Bz) + * - Electric field: 3 components (Ex, Ey, Ez) at indices 3-5 + * - Electromagnetic field: 6 components (Bx, By, Bz, Ex, Ey, Ez) + */ +class PyG4Field : public G4Field { +public: + // Inherit the constructors + using G4Field::G4Field; + + void GetFieldValue(const G4double Point[4], + G4double *fieldArr) const override { + py::gil_scoped_acquire gil; + // Always initialize output to a safe value. + for (size_t i = 0; i < G4Field::MAX_NUMBER_OF_COMPONENTS; ++i) { + fieldArr[i] = 0.; + } + + // Convert Point to Python list + py::list pyPoint; + pyPoint.append(Point[0]); + pyPoint.append(Point[1]); + pyPoint.append(Point[2]); + pyPoint.append(Point[3]); + + // Try to get the Python override + py::function override = py::get_override(this, "GetFieldValue"); + if (override) { + // Call Python implementation and expect a list back + py::object result = override(pyPoint); + if (!result.is_none()) { + py::list field_list = result.cast(); + size_t n = py::len(field_list); + for (size_t i = 0; i < n && i < G4Field::MAX_NUMBER_OF_COMPONENTS; + ++i) { + fieldArr[i] = field_list[i].cast(); + } + } + } + } + + G4bool DoesFieldChangeEnergy() const override { + PYBIND11_OVERRIDE_PURE(G4bool, G4Field, DoesFieldChangeEnergy); + } +}; + +void init_G4Field(py::module &m) { + + py::class_>( + m, "G4Field") + + .def(py::init()) + .def(py::init<>()) + + .def("DoesFieldChangeEnergy", &G4Field::DoesFieldChangeEnergy) + .def("IsGravityActive", &G4Field::IsGravityActive) + .def("SetGravityActive", &G4Field::SetGravityActive) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4FieldManager.cpp b/core/opengate_core/g4_bindings/pyG4FieldManager.cpp new file mode 100644 index 000000000..5ec7d2292 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4FieldManager.cpp @@ -0,0 +1,66 @@ +/* -------------------------------------------------- + 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 "G4ChordFinder.hh" +#include "G4Field.hh" +#include "G4FieldManager.hh" +#include "G4MagneticField.hh" +#include "G4Track.hh" + +void init_G4FieldManager(py::module &m) { + py::class_>( + m, "G4FieldManager") + + .def(py::init()) + .def(py::init()) + + .def("SetDetectorField", &G4FieldManager::SetDetectorField) + .def("ProposeDetectorField", &G4FieldManager::ProposeDetectorField) + .def("ChangeDetectorField", &G4FieldManager::ChangeDetectorField) + .def("GetDetectorField", &G4FieldManager::GetDetectorField, + py::return_value_policy::reference_internal) + + .def("DoesFieldExist", &G4FieldManager::DoesFieldExist) + + .def("CreateChordFinder", &G4FieldManager::CreateChordFinder) + .def("SetChordFinder", &G4FieldManager::SetChordFinder) + .def("GetChordFinder", + py::overload_cast<>(&G4FieldManager::GetChordFinder), + py::return_value_policy::reference_internal) + + .def("ConfigureForTrack", &G4FieldManager::ConfigureForTrack) + + // Cannot be compiled on windows because of the static member + // G4FieldManager::GlobalFieldManager + //.def("SetGlobalFieldManager", &G4FieldManager::SetGlobalFieldManager) + //.def("GetGlobalFieldManager", &G4FieldManager::GetGlobalFieldManager, + // py::return_value_policy::reference_internal) + + .def("GetDeltaIntersection", &G4FieldManager::GetDeltaIntersection) + .def("GetDeltaOneStep", &G4FieldManager::GetDeltaOneStep) + .def("SetAccuraciesWithDeltaOneStep", + &G4FieldManager::SetAccuraciesWithDeltaOneStep) + .def("SetDeltaOneStep", &G4FieldManager::SetDeltaOneStep) + .def("SetDeltaIntersection", &G4FieldManager::SetDeltaIntersection) + + .def("GetMinimumEpsilonStep", &G4FieldManager::GetMinimumEpsilonStep) + .def("SetMinimumEpsilonStep", &G4FieldManager::SetMinimumEpsilonStep) + .def("GetMaximumEpsilonStep", &G4FieldManager::GetMaximumEpsilonStep) + .def("SetMaximumEpsilonStep", &G4FieldManager::SetMaximumEpsilonStep) + + .def("DoesFieldChangeEnergy", &G4FieldManager::DoesFieldChangeEnergy) + .def("SetFieldChangesEnergy", &G4FieldManager::SetFieldChangesEnergy) + + .def("GetMaxAcceptedEpsilon", &G4FieldManager::GetMaxAcceptedEpsilon) + .def("SetMaxAcceptedEpsilon", &G4FieldManager::SetMaxAcceptedEpsilon) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4MagErrorStepper.cpp b/core/opengate_core/g4_bindings/pyG4MagErrorStepper.cpp new file mode 100644 index 000000000..4a90d5928 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4MagErrorStepper.cpp @@ -0,0 +1,22 @@ +/* -------------------------------------------------- + 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 "G4MagErrorStepper.hh" +#include "G4MagIntegratorStepper.hh" + +void init_G4MagErrorStepper(py::module &m) { + + py::class_>( + m, "G4MagErrorStepper") + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4MagInt_Driver.cpp b/core/opengate_core/g4_bindings/pyG4MagInt_Driver.cpp new file mode 100644 index 000000000..9745d57ff --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4MagInt_Driver.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 + +namespace py = pybind11; + +#include "G4MagIntegratorDriver.hh" +#include "G4MagIntegratorStepper.hh" +#include "G4VIntegrationDriver.hh" + +void init_G4MagInt_Driver(py::module &m) { + + py::class_>(m, + "G4MagInt_Driver") + + .def(py::init()) + + .def("GetHmin", &G4MagInt_Driver::GetHmin) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4MagIntegratorStepper.cpp b/core/opengate_core/g4_bindings/pyG4MagIntegratorStepper.cpp new file mode 100644 index 000000000..a71895784 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4MagIntegratorStepper.cpp @@ -0,0 +1,28 @@ +/* -------------------------------------------------- + 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 "G4EquationOfMotion.hh" +#include "G4MagIntegratorStepper.hh" + +void init_G4MagIntegratorStepper(py::module &m) { + + py::class_>( + m, "G4MagIntegratorStepper") + + .def("GetNumberOfVariables", + &G4MagIntegratorStepper::GetNumberOfVariables) + .def("GetNumberOfStateVariables", + &G4MagIntegratorStepper::GetNumberOfStateVariables) + .def("IntegratorOrder", &G4MagIntegratorStepper::IntegratorOrder) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4Mag_EqRhs.cpp b/core/opengate_core/g4_bindings/pyG4Mag_EqRhs.cpp new file mode 100644 index 000000000..a8ac8faf1 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4Mag_EqRhs.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 + +namespace py = pybind11; + +#include "G4EquationOfMotion.hh" +#include "G4Mag_EqRhs.hh" + +void init_G4Mag_EqRhs(py::module &m) { + + py::class_>(m, "G4Mag_EqRhs") + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4Mag_UsualEqRhs.cpp b/core/opengate_core/g4_bindings/pyG4Mag_UsualEqRhs.cpp new file mode 100644 index 000000000..6440095b3 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4Mag_UsualEqRhs.cpp @@ -0,0 +1,29 @@ +/* -------------------------------------------------- + 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 "G4EquationOfMotion.hh" +#include "G4Mag_EqRhs.hh" +#include "G4Mag_UsualEqRhs.hh" +#include "G4MagneticField.hh" + +void init_G4Mag_UsualEqRhs(py::module &m) { + + py::class_>( + m, "G4Mag_UsualEqRhs") + + .def(py::init()) + + .def("EvaluateRhsGivenB", &G4Mag_UsualEqRhs::EvaluateRhsGivenB) + .def("SetChargeMomentumMass", &G4Mag_UsualEqRhs::SetChargeMomentumMass) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4MagneticField.cpp b/core/opengate_core/g4_bindings/pyG4MagneticField.cpp new file mode 100644 index 000000000..9257d42d5 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4MagneticField.cpp @@ -0,0 +1,75 @@ +/* -------------------------------------------------- + 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 "G4Field.hh" +#include "G4MagneticField.hh" + +/* + * Trampoline class to allow Python to override GetFieldValue for magnetic + * fields. Inherits from G4MagneticField which already implements + * DoesFieldChangeEnergy() = false. + * + * The GetFieldValue method receives a point (x, y, z, t) and must return + * the magnetic field vector (Bx, By, Bz) at that point. + */ +class PyG4MagneticField : public G4MagneticField { +public: + using G4MagneticField::G4MagneticField; + + void GetFieldValue(const G4double Point[4], G4double *Bfield) const override { + // Always initialize output to a safe value. + Bfield[0] = 0.; + Bfield[1] = 0.; + Bfield[2] = 0.; + + py::gil_scoped_acquire gil; + + // Convert Point to Python list + py::list pyPoint; + pyPoint.append(Point[0]); + pyPoint.append(Point[1]); + pyPoint.append(Point[2]); + pyPoint.append(Point[3]); + + // Try to get the Python override + py::function override = py::get_override(this, "GetFieldValue"); + if (override) { + py::object result = override(pyPoint); + if (!result.is_none()) { + py::sequence bfield_seq = result.cast(); + + if (py::len(bfield_seq) != 3) { + throw std::invalid_argument( + "GetFieldValue for G4MagneticField must return exactly 3 " + "components [Bx, By, Bz]"); + } + + Bfield[0] = bfield_seq[0].cast(); + Bfield[1] = bfield_seq[1].cast(); + Bfield[2] = bfield_seq[2].cast(); + } + } + } +}; + +void init_G4MagneticField(py::module &m) { + + py::class_>(m, + "G4MagneticField") + + .def(py::init<>()) + + .def("DoesFieldChangeEnergy", &G4MagneticField::DoesFieldChangeEnergy) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4QuadrupoleMagField.cpp b/core/opengate_core/g4_bindings/pyG4QuadrupoleMagField.cpp new file mode 100644 index 000000000..d59f0ce95 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4QuadrupoleMagField.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 "G4QuadrupoleMagField.hh" + +void init_G4QuadrupoleMagField(py::module &m) { + + py::class_>( + m, "G4QuadrupoleMagField") + + .def(py::init()) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4UniformElectricField.cpp b/core/opengate_core/g4_bindings/pyG4UniformElectricField.cpp new file mode 100644 index 000000000..1ea06e1e0 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4UniformElectricField.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 "G4ElectricField.hh" +#include "G4UniformElectricField.hh" + +void init_G4UniformElectricField(py::module &m) { + + py::class_>( + m, "G4UniformElectricField") + + .def(py::init()) + .def(py::init()) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4UniformMagField.cpp b/core/opengate_core/g4_bindings/pyG4UniformMagField.cpp new file mode 100644 index 000000000..a384dcf78 --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4UniformMagField.cpp @@ -0,0 +1,28 @@ +/* -------------------------------------------------- + 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 "G4UniformMagField.hh" + +void init_G4UniformMagField(py::module &m) { + + py::class_>( + m, "G4UniformMagField") + + .def(py::init()) + + .def("SetFieldValue", py::overload_cast( + &G4UniformMagField::SetFieldValue)) + .def("GetConstantFieldValue", &G4UniformMagField::GetConstantFieldValue) + + ; +} diff --git a/core/opengate_core/g4_bindings/pyG4VIntegrationDriver.cpp b/core/opengate_core/g4_bindings/pyG4VIntegrationDriver.cpp new file mode 100644 index 000000000..5bf9c613a --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4VIntegrationDriver.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 + +namespace py = pybind11; + +#include "G4VIntegrationDriver.hh" + +void init_G4VIntegrationDriver(py::module &m) { + + py::class_>( + m, "G4VIntegrationDriver") + + .def("SetEquationOfMotion", &G4VIntegrationDriver::SetEquationOfMotion) + .def("GetEquationOfMotion", &G4VIntegrationDriver::GetEquationOfMotion, + py::return_value_policy::reference_internal) + .def("SetVerboseLevel", &G4VIntegrationDriver::SetVerboseLevel) + .def("GetVerboseLevel", &G4VIntegrationDriver::GetVerboseLevel) + + ; +} diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index 3737dea8f..9252ffd6f 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -173,6 +173,42 @@ void init_G4Region(py::module &); void init_G4RegionStore(py::module &); +void init_G4FieldManager(py::module &); + +void init_G4Field(py::module &); + +void init_G4MagneticField(py::module &); + +void init_G4ElectroMagneticField(py::module &); + +void init_G4ElectricField(py::module &); + +void init_G4UniformMagField(py::module &); + +void init_G4QuadrupoleMagField(py::module &); + +void init_G4UniformElectricField(py::module &); + +void init_G4EquationOfMotion(py::module &); + +void init_G4Mag_EqRhs(py::module &); + +void init_G4Mag_UsualEqRhs(py::module &); + +void init_G4EqMagElectricField(py::module &); + +void init_G4MagIntegratorStepper(py::module &); + +void init_G4MagErrorStepper(py::module &); + +void init_G4ClassicalRK4(py::module &); + +void init_G4VIntegrationDriver(py::module &); + +void init_G4MagInt_Driver(py::module &); + +void init_G4ChordFinder(py::module &); + // geometry/solids void init_G4Box(py::module &); @@ -502,6 +538,25 @@ PYBIND11_MODULE(opengate_core, m) { init_G4Region(m); init_G4RegionStore(m); + init_G4FieldManager(m); + init_G4Field(m); + init_G4MagneticField(m); + init_G4ElectroMagneticField(m); + init_G4ElectricField(m); + init_G4UniformMagField(m); + init_G4QuadrupoleMagField(m); + init_G4UniformElectricField(m); + init_G4EquationOfMotion(m); + init_G4Mag_EqRhs(m); + init_G4Mag_UsualEqRhs(m); + init_G4EqMagElectricField(m); + init_G4MagIntegratorStepper(m); + init_G4MagErrorStepper(m); + init_G4ClassicalRK4(m); + init_G4VIntegrationDriver(m); + init_G4MagInt_Driver(m); + init_G4ChordFinder(m); + init_G4Box(m); init_G4Ellipsoid(m); init_G4Polyhedra(m); diff --git a/docs/source/user_guide/index.rst b/docs/source/user_guide/index.rst index 10c7c099f..55d7eda44 100644 --- a/docs/source/user_guide/index.rst +++ b/docs/source/user_guide/index.rst @@ -19,6 +19,7 @@ User Guide user_guide_physics user_guide_sources user_guide_actors + user_guide_fields .. toctree:: :caption: Detailed description of components @@ -28,6 +29,7 @@ User Guide user_guide_reference_volumes user_guide_reference_actors user_guide_reference_sources + user_guide_reference_fields user_guide_reference_filters user_guide_advanced user_guide_optical_photons diff --git a/docs/source/user_guide/user_guide_fields.rst b/docs/source/user_guide/user_guide_fields.rst new file mode 100644 index 000000000..247bd0d06 --- /dev/null +++ b/docs/source/user_guide/user_guide_fields.rst @@ -0,0 +1,272 @@ +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. + +Quick example +------------- + +.. code-block:: python + + import opengate as gate + from opengate.geometry import fields + + sim = gate.Simulation() + tesla = gate.g4_units.tesla + cm = gate.g4_units.cm + + box = sim.add_volume("Box", "magnet") + box.size = [50 * cm, 50 * cm, 50 * cm] + box.material = "G4_Galactic" + + field = fields.UniformMagneticField(name="B_field_for_box") + field.field_vector = [0, 1 * tesla, 0] + box.add_field(field) + + +Available field types +--------------------- + +Magnetic fields +~~~~~~~~~~~~~~~ + +**UniformMagneticField** -- Constant magnetic field throughout the volume. Uses the native Geant4 ``G4UniformMagField`` under the hood. + +.. code-block:: python + + field = fields.UniformMagneticField(name="B_uniform") + field.field_vector = [0, 0, 2 * tesla] # [Bx, By, Bz] + box.add_field(field) + +**QuadrupoleMagneticField** -- Quadrupole magnetic field defined via a field gradient. Uses the native Geant4 ``G4QuadrupoleMagField`` under the hood. + +.. code-block:: python + + field = fields.QuadrupoleMagneticField(name="B_quad") + 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): + # Spatially varying field + return [0, (1 + x * z / m**2) * tesla, 0] + + field = fields.CustomMagneticField(name="B_custom") + field.field_function = my_B_field + box.add_field(field) + +.. warning:: Performance warning + + 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`. + + +Electric fields +~~~~~~~~~~~~~~~ + +**UniformElectricField** -- Constant electric field. Uses the native Geant4 ``G4UniformElectricField`` under the hood. + +.. code-block:: python + + volt = gate.g4_units.volt + m = gate.g4_units.m + + field = fields.UniformElectricField(name="E_uniform") + field.field_vector = [1e6 * volt / m, 0, 0] # [Ex, Ey, Ez] + 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] + + field = fields.CustomElectricField(name="E_custom") + field.field_function = my_E_field + box.add_field(field) + + +Electromagnetic fields +~~~~~~~~~~~~~~~~~~~~~~ + +Combined magnetic and electric fields. + +**CustomElectroMagneticField** -- Arbitrary combined field. The callback must return all six components ``[Bx, By, Bz, Ex, Ey, Ez]``. + +.. code-block:: python + + def my_EM_field(x, y, z, t): + return [0, 1 * tesla, 0, 1e6 * volt / m, 0, 0] + + field = fields.CustomElectroMagneticField(name="EM_custom") + field.field_function = my_EM_field + box.add_field(field) + + +Integration 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. + +.. list-table:: + :header-rows: 1 + :widths: 30 15 55 + + * - Parameter + - Default + - Description + * - ``step_minimum`` + - 0.01 mm + - Minimum step size for the chord finder. + * - ``delta_chord`` + - 0.001 mm + - Maximum sagitta (miss distance between the chord approximation and the true curved trajectory). + * - ``delta_one_step`` + - 0.001 mm + - Positional accuracy per integration step. + * - ``delta_intersection`` + - 0.0001 mm + - Positional accuracy at volume boundaries. + * - ``min_epsilon_step`` + - 1e-7 + - Minimum relative integration accuracy. + * - ``max_epsilon_step`` + - 1e-5 + - Maximum relative integration accuracy. + +Example: + +.. code-block:: python + + field = fields.UniformMagneticField(name="B") + field.field_vector = [0, 0, 1 * tesla] + field.delta_chord = 0.01 * mm + field.step_minimum = 0.1 * mm + box.add_field(field) + + +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. + +.. code-block:: python + + box = sim.add_volume("Box", "my_box") + # ...configure box and define field + box.add_field(field) + + +**One field per volume.** Attempting to attach a second field raises an error. + +**Shared fields.** The same field object can be attached to multiple volumes. This is useful when several volumes should have the same field configuration: + +.. code-block:: python + + field = fields.UniformMagneticField(name="B_shared") + field.field_vector = [0, 0, 1 * tesla] + box1.add_field(field) + box2.add_field(field) + +**Unique names.** Field names must be unique across the simulation. Two different field objects with the same name will raise an error. + +**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. + +.. code-block:: python + + outer = sim.add_volume("Box", "outer") + inner = sim.add_volume("Box", "inner") + inner.mother = "outer" + + field_outer = fields.UniformMagneticField(name="B_outer") + field_outer.field_vector = [0, 0, 1 * tesla] + outer.add_field(field_outer) + # "inner" inherits B_outer automatically. + + # To override it inside the daughter, attach a different field: + field_inner = fields.UniformMagneticField(name="B_inner") + field_inner.field_vector = [0, 2 * tesla, 0] + inner.add_field(field_inner) + # Now "inner" uses B_inner; "outer" still uses B_outer outside of "inner". + + +Visualizing fields +------------------ + +Geant4 can overlay field arrows on the visualization of the geometry. The Geant4 command is ``/vis/scene/add/magneticField `` (or ``/vis/scene/add/electricField`` for electric fields). It can be passed to GATE via the ``visu_commands`` parameter: + +.. code-block:: python + + sim.visu = True + sim.visu_type = "qt" + sim.visu_commands.append("/vis/scene/add/magneticField 20 lightArrow") + + +.. _user_guide_fields_performance: + +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. + +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. + +.. 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 +.. -------------- + +.. Non-custom fields are fully serializable via GATE's ``to_dictionary()`` / ``from_dictionary()`` mechanism: + +.. .. code-block:: python + +.. d = sim.to_dictionary() +.. sim2 = gate.Simulation() +.. sim2.from_dictionary(d) +.. # sim2 now has the same fields attached to the same volumes + +.. Custom fields (those with a ``field_function`` callback) **cannot** be serialized. + + +Examples +-------- + +The field implementation is covered by the ``test099_fields_*`` tests in ``opengate/tests/src/geometry/``. These tests can be used as examples of how to define and use fields in GATE. They include: + +- ``test099_fields_analytical_B`` -- Uniform B field vs analytical cyclotron radius. +- ``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_serialization`` -- Round-trip serialization for all non-custom types. +- ``test099_fields_api`` -- API guards. + + +Class reference +---------------- + +.. autoclass:: opengate.geometry.fields.FieldBase + :members: + +.. autoclass:: opengate.geometry.fields.UniformMagneticField + :members: + +.. autoclass:: opengate.geometry.fields.QuadrupoleMagneticField + :members: + +.. autoclass:: opengate.geometry.fields.CustomMagneticField + :members: + +.. autoclass:: opengate.geometry.fields.UniformElectricField + :members: + +.. autoclass:: opengate.geometry.fields.CustomElectricField + :members: + +.. autoclass:: opengate.geometry.fields.CustomElectroMagneticField + :members: diff --git a/docs/source/user_guide/user_guide_reference_fields.rst b/docs/source/user_guide/user_guide_reference_fields.rst new file mode 100644 index 000000000..9aedbfeae --- /dev/null +++ b/docs/source/user_guide/user_guide_reference_fields.rst @@ -0,0 +1,25 @@ +.. _fields-reference-label: + +Details: Electromagnetic fields +================================ + +[TO BE COMPLETED] + +Base class reference +-------------------- + +.. autoclass:: opengate.geometry.fields.FieldBase + :members: + :show-inheritance: + +.. autoclass:: opengate.geometry.fields.MagneticField + :members: + :show-inheritance: + +.. autoclass:: opengate.geometry.fields.ElectroMagneticField + :members: + :show-inheritance: + +.. autoclass:: opengate.geometry.fields.ElectricField + :members: + :show-inheritance: diff --git a/opengate/__init__.py b/opengate/__init__.py index 8fbe32c0b..8ced0adbf 100644 --- a/opengate/__init__.py +++ b/opengate/__init__.py @@ -135,6 +135,7 @@ def check_ld_preload(e): import opengate.geometry.solids import opengate.geometry.utility import opengate.geometry.volumes +import opengate.geometry.fields import opengate.actors import opengate.contrib diff --git a/opengate/engines.py b/opengate/engines.py index f4a5eccf3..365153e12 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -861,6 +861,14 @@ def ConstructSDandField(self): self.volume_manager.world_volume.name, ) + 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) + volume_obj.g4_field_manager = field.create_field_manager() + volume_obj.g4_logical_volume.SetFieldManager( + volume_obj.g4_field_manager, True + ) + def get_volume(self, name): return self.volume_manager.get_volume(name) diff --git a/opengate/geometry/fields.py b/opengate/geometry/fields.py new file mode 100644 index 000000000..a487bcb91 --- /dev/null +++ b/opengate/geometry/fields.py @@ -0,0 +1,553 @@ +import opengate_core as g4 + +from ..base import GateObject, process_cls +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) +# ! - 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. +# ! ============================= + + +class FieldBase(GateObject): + """Base class for electric and magnetic fields.""" + + # hints for IDE + step_minimum: float + delta_chord: float + delta_one_step: float + delta_intersection: float + min_epsilon_step: float + max_epsilon_step: float + + user_info_defaults = { + "step_minimum": ( + 1e-2 * g4_units.mm, + { + "doc": "Minimum step size for the chord finder.", + }, + ), + "delta_chord": ( + 1e-3 * g4_units.mm, + { + "doc": "Maximum miss distance between chord and curved trajectory.", + }, + ), + "delta_one_step": ( + 1e-3 * g4_units.mm, + { + "doc": "Positional accuracy per integration step.", + }, + ), + "delta_intersection": ( + 1e-4 * g4_units.mm, + { + "doc": "Positional accuracy at volume boundaries.", + }, + ), + "min_epsilon_step": ( + 1e-7, + { + "doc": "Minimum relative integration accuracy.", + }, + ), + "max_epsilon_step": ( + 1e-5, + { + "doc": "Maximum relative integration accuracy.", + }, + ), + } + + def __init__(self, *args, **kwargs) -> None: + super().__init__(*args, **kwargs) + + self.g4_field = None + self._g4_runtime_objects = [] + + self.attached_to = [] + self._field_changes_energy = False + + @property + def field_changes_energy(self) -> bool: + """Whether the field changes particle energy (False for magnetic, True for others).""" + return self._field_changes_energy + + 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) + # 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_field_manager(self) -> g4.G4FieldManager: + """Construct the field and return a configured G4FieldManager.""" + # Create the G4 field object (subclass responsibility) + self._create_field() + + # Create equation of motion, stepper, driver, chord finder TODO: allow user choice + 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) + ) + 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) + 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 + ) + 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, + } + ) + + 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() + + +class UniformMagneticField(MagneticField): + """Uniform magnetic field with constant field vector.""" + + # hints for IDE + field_vector: list + + user_info_defaults = { + "field_vector": ( + [0, 0, 0], + { + "doc": "Field vector [Bx, By, Bz]. Each component in magnetic field strength units (e.g., Tesla).", + }, + ), + } + + 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)) + + +class QuadrupoleMagneticField(MagneticField): + """Quadrupole 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 __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) + + +class CustomMagneticField(MagneticField): + """Custom magnetic field defined by a Python callback function.""" + + # hints for IDE + field_function: callable + + user_info_defaults = { + "field_function": ( + None, + { + "doc": "Python function that takes [x, y, z, t] and returns [Bx, By, Bz].", + }, + ), + } + + def __init__(self, *args, **kwargs) -> None: + super().__init__(*args, **kwargs) + + def _create_field(self) -> None: + """Create the custom magnetic field using the Python trampoline.""" + 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) + 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__() + inner_self._callback = callback + + def GetFieldValue(inner_self, point): + return inner_self._callback(*point) + + self.g4_field = _PyMagneticField(self.field_function) + + def to_dictionary(self): + raise NotImplementedError( + "Custom fields with Python callbacks cannot be serialized." + ) + + +class ElectroMagneticField(FieldBase): + """Base class for electromagnetic fields.""" + + 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_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, + } + ) + + 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() + + +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 + + +class UniformElectricField(ElectricField): + """Uniform electric field with constant field vector.""" + + # hints for IDE + field_vector: list + + user_info_defaults = { + "field_vector": ( + [0, 0, 0], + { + "doc": "Field vector [Ex, Ey, Ez] in Geant4 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)) + + +class CustomElectricField(ElectricField): + """Custom electric field defined by a Python callback function.""" + + # hints for IDE + field_function: callable + + user_info_defaults = { + "field_function": ( + None, + { + "doc": "Python function that takes [x, y, z, t] and returns [Ex, Ey, Ez].", + }, + ), + } + + 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]" + ) + + # Create a custom G4ElectricField subclass that calls our Python function + class _PyElectricField(g4.G4ElectricField): + def __init__(inner_self, callback): + super().__init__() + inner_self._callback = callback + + def GetFieldValue(inner_self, point): + return inner_self._callback(*point) + + def DoesFieldChangeEnergy(inner_self): + return True + + self.g4_field = _PyElectricField(self.field_function) + + def to_dictionary(self): + raise NotImplementedError( + "Custom fields with Python callbacks cannot be serialized." + ) + + +# 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.""" + + # hints for IDE + magnetic_field_vector: list + electric_field_vector: list + + user_info_defaults = { + "magnetic_field_vector": ( + [0, 0, 0], + { + "doc": "Magnetic field vector [Bx, By, Bz] in Geant4 units.", + }, + ), + "electric_field_vector": ( + [0, 0, 0], + { + "doc": "Electric field vector [Ex, Ey, Ez] in Geant4 units.", + }, + ), + } + + 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]) + + +class CustomElectroMagneticField(ElectroMagneticField): + """Custom electromagnetic field defined by a Python callback function.""" + + # hints for IDE + field_function: callable + + user_info_defaults = { + "field_function": ( + None, + { + "doc": "Python function that takes [x, y, z, t] and returns [Bx, By, Bz, Ex, Ey, Ez].", + }, + ), + } + + 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]" + ) + + # Create a custom G4ElectroMagneticField subclass that calls our Python function + class _PyElectroMagneticField(g4.G4ElectroMagneticField): + def __init__(inner_self, callback): + super().__init__() + inner_self._callback = callback + + def GetFieldValue(inner_self, point): + return inner_self._callback(*point) + + def DoesFieldChangeEnergy(inner_self): + return True + + self.g4_field = _PyElectroMagneticField(self.field_function) + + def to_dictionary(self): + raise NotImplementedError( + "Custom fields with Python callbacks cannot be serialized." + ) + + +field_types = { + "UniformMagneticField": UniformMagneticField, + "QuadrupoleMagneticField": QuadrupoleMagneticField, + "CustomMagneticField": CustomMagneticField, + "UniformElectricField": UniformElectricField, + "CustomElectricField": CustomElectricField, + "UniformElectroMagneticField": UniformElectroMagneticField, + "CustomElectroMagneticField": CustomElectroMagneticField, +} + +process_cls(FieldBase) +process_cls(MagneticField) +process_cls(UniformMagneticField) +process_cls(QuadrupoleMagneticField) +process_cls(CustomMagneticField) +process_cls(ElectroMagneticField) +process_cls(ElectricField) +process_cls(UniformElectricField) +process_cls(CustomElectricField) +process_cls(UniformElectroMagneticField) +process_cls(CustomElectroMagneticField) diff --git a/opengate/geometry/volumes.py b/opengate/geometry/volumes.py index 5df29d539..524336bbb 100644 --- a/opengate/geometry/volumes.py +++ b/opengate/geometry/volumes.py @@ -23,6 +23,7 @@ ) from ..decorators import requires_fatal, requires_attribute_fatal from ..definitions import __world_name__, __gate_list_objects__ +from .fields import FieldBase from ..actors.dynamicactors import ( VolumeImageChanger, VolumeTranslationChanger, @@ -185,6 +186,13 @@ class VolumeBase(DynamicGateObject, NodeMixin): "type": bool, }, ), + "field": ( + None, + { + "doc": "Name of the field attached to this volume.", + "type": str, + }, + ), } def __init__(self, *args, template=None, **kwargs): @@ -216,6 +224,10 @@ def __init__(self, *args, template=None, **kwargs): # this list contains all physical volumes (in case of repeated volume) self.g4_physical_volumes = [] self.g4_material = None + self.g4_field_manager = None + + # Field attached to this volume (only one allowed) + self.field = None def close(self): self.release_g4_references() @@ -228,6 +240,7 @@ def release_g4_references(self): self.g4_vis_attributes = None self.g4_physical_volumes = [] self.g4_material = None + self.g4_field_manager = None def __getstate__(self): return_dict = super().__getstate__() @@ -237,6 +250,7 @@ def __getstate__(self): return_dict["g4_vis_attributes"] = None return_dict["g4_physical_volumes"] = [] return_dict["g4_material"] = None + return_dict["g4_field_manager"] = None return_dict["volume_engine"] = None return_dict["_is_constructed"] = False return return_dict @@ -531,6 +545,23 @@ def set_min_range(self, min_range): self.name, min_range ) + @requires_fatal("volume_manager") + def add_field(self, field: FieldBase): + if self.field is not None: + fatal( + f"Volume '{self.name}' already has a field attached ('{self.field}'). " + f"A volume can only have one field. Remove the existing field first." + ) + existing = self.volume_manager.fields.get(field.name) + if existing is not None and existing is not field: + fatal( + f"A field named '{field.name}' is already registered in the volume manager. " + f"Field names must be unique, please choose a different name for this field." + ) + self.field = field.name + field.attached_to.append(self.name) + self.volume_manager.fields[field.name] = field + class RepeatableVolume(VolumeBase): def get_repetition_name_from_index(self, index): diff --git a/opengate/managers.py b/opengate/managers.py index 3b7809949..206e84738 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -18,6 +18,7 @@ from .engines import SimulationEngine from .exception import fatal, warning, GateDeprecationError, GateImplementationError from .geometry.materials import MaterialDatabase +from .geometry.fields import FieldBase, field_types from .utility import ( g4_units, @@ -1159,6 +1160,9 @@ def __init__(self, simulation, *args, **kwargs) -> None: # Store them to init them later self.solid_with_texture_init = [] + # fields + self.fields: dict[str, FieldBase] = {} + def reset(self): self.__init__(self.simulation) @@ -1176,11 +1180,23 @@ def to_dictionary(self): d = super().to_dictionary() d["volumes"] = dict([(k, v.to_dictionary()) for k, v in self.volumes.items()]) d["parallel_world_volumes"] = list(self.parallel_world_volumes.keys()) + d["fields"] = dict([(k, v.to_dictionary()) for k, v in self.fields.items()]) return d def from_dictionary(self, d): self.reset() super().from_dictionary(d) + # Restore fields before volumes so that volume references are valid + for k, v in d.get("fields", {}).items(): + field_type = v["object_type"] + if field_type not in field_types: + fatal( + f"Unknown field type '{field_type}'. " + f"Known types are: {list(field_types.keys())}." + ) + field = field_types[field_type](name=v["user_info"]["name"]) + field.from_dictionary(v) + self.fields[field.name] = field # First create all volumes for k, v in d["volumes"].items(): # the world volume is always created in __init__ diff --git a/opengate/tests/src/geometry/test099_fields_analytical_B.py b/opengate/tests/src/geometry/test099_fields_analytical_B.py new file mode 100644 index 000000000..39415c484 --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_analytical_B.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python3 +""" +Test 099 — Uniform magnetic field: analytical validation. + +A proton enters a region with uniform B along Y. +Checks: + 1. Circular trajectory in XZ plane (cyclotron radius) + 2. No deflection in Y + 3. Energy conservation (B does no work) +""" + +import numpy as np +import uproot + +from opengate.geometry import fields +from opengate.tests import utility + +from test099_fields_helpers import ( + g4_tesla, + g4_mm, + g4_MeV, + g4_eplus, + PROTON_MASS, + cyclotron_radius, + build_field_simulation, +) + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, output_folder="test099_fields") + + By = 5 * g4_tesla + T = 10 * g4_MeV + + field = fields.UniformMagneticField(name="B_uniform") + field.field_vector = [0, By, 0] + + sim, phsp = build_field_simulation( + field, + kinetic_energy=T, + phsp_output_filename=paths.output / "test099_analytical_B.root", + output_dir=paths.output, + ) + sim.run() + + # Read results + df = uproot.open(str(paths.output / "test099_analytical_B.root"))["phsp;1"].arrays( + library="pd" + ) + + # Analytical expectation + r = cyclotron_radius(T, By, PROTON_MASS, 1 * g4_eplus) + + # Entry point & circle centre + box_half_z = 50 * g4_mm / 2 * 10 + x0, y0, z0 = 0.0, 0.0, -box_half_z + cx, cz = x0 - r, z0 + + x_exit = df["PostPosition_X"].values + y_exit = df["PostPosition_Y"].values + z_exit = df["PostPosition_Z"].values + KE_exit = df["PostKineticEnergy"].values + + # Check 1: circular trajectory + r_TOL = 1e-3 * g4_mm + residual = np.sqrt((x_exit - cx) ** 2 + (z_exit - cz) ** 2) - r + is_ok_circular = np.all(np.abs(residual) < r_TOL) + print(f"Circle-fit residual: {residual}") + print(f"Circular trajectory OK: {is_ok_circular}") + + # Check 2: no Y deflection + is_ok_y = np.all(np.abs(y_exit - y0) < r_TOL) + print(f"No Y deflection OK: {is_ok_y}") + + # Check 3: energy conservation + e_TOL = 0.01 * g4_MeV + is_ok_energy = np.all(np.abs(KE_exit - T) < e_TOL) + print(f"Energy conservation OK: {is_ok_energy} (dKE = {KE_exit - T})") + + is_ok = is_ok_circular and is_ok_y and is_ok_energy + utility.test_ok(is_ok) diff --git a/opengate/tests/src/geometry/test099_fields_analytical_E.py b/opengate/tests/src/geometry/test099_fields_analytical_E.py new file mode 100644 index 000000000..9203b5d28 --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_analytical_E.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python3 +""" +Test 099 — Uniform electric field: analytical validation. + +A proton enters a region with uniform E along X. +Checks: + 1. Displacement consistent with energy gain: dx = dKE / (qE) + 2. No deflection in Y + 3. Energy gain matches work done: dKE = q * E * dx +""" + +import numpy as np +import uproot + +from opengate.geometry import fields +from opengate.tests import utility + +from test099_fields_helpers import ( + g4_m, + g4_mm, + g4_MeV, + g4_volt, + g4_eplus, + PROTON_MASS, + build_field_simulation, +) + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, output_folder="test099_fields") + + Ex = 1e8 * g4_volt / g4_m + T = 10 * g4_MeV + + field = fields.UniformElectricField(name="E_uniform") + field.field_vector = [Ex, 0, 0] + + sim, phsp = build_field_simulation( + field, + kinetic_energy=T, + phsp_output_filename=paths.output / "test099_analytical_E.root", + output_dir=paths.output, + ) + sim.run() + + df = uproot.open(str(paths.output / "test099_analytical_E.root"))["phsp;1"].arrays( + library="pd" + ) + + q = 1 * g4_eplus + x_exit = df["PostPosition_X"].values + y_exit = df["PostPosition_Y"].values + KE_exit = df["PostKineticEnergy"].values + r_TOL = 0.01 * g4_mm + e_TOL = 0.01 * g4_MeV + + # Check 1: x displacement consistent with energy gain + x_from_energy = (KE_exit - T) / (q * Ex) + is_ok_x = np.all(np.abs(x_exit - x_from_energy) < r_TOL) + print(f"x consistency OK: {is_ok_x} (residual = {x_exit - x_from_energy})") + + # Check 2: no Y deflection + is_ok_y = np.all(np.abs(y_exit) < r_TOL) + print(f"No Y deflection OK: {is_ok_y}") + + # Check 3: energy gain = qE * dx + KE_expected = T + q * Ex * x_exit / g4_MeV + is_ok_energy = np.all(np.abs(KE_exit - KE_expected) < e_TOL) + print(f"Energy gain OK: {is_ok_energy} (residual = {KE_exit - KE_expected})") + + is_ok = is_ok_x and is_ok_y and is_ok_energy + utility.test_ok(is_ok) diff --git a/opengate/tests/src/geometry/test099_fields_api.py b/opengate/tests/src/geometry/test099_fields_api.py new file mode 100644 index 000000000..294b77e1a --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_api.py @@ -0,0 +1,149 @@ +#!/usr/bin/env python3 +""" +Test 099 — Field API tests. + +Checks the following behaviours: + 1. add_field() before volume is added to simulation -> fatal + 2. Duplicate field on same volume -> fatal + 3. Two different fields with same name -> fatal + 4. Same field on two volumes -> OK + 5. Custom field with non-callable -> TypeError + 6. Custom field with wrong component count -> ValueError +""" + +import opengate as gate +from opengate.geometry import fields +from opengate.geometry.volumes import VolumeBase +from opengate.tests import utility + +from test099_fields_helpers import g4_tesla, g4_m, g4_cm + + +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}") + return False + except Exception as e: + print(f" OK — got {type(e).__name__}: {error_msg}") + return True + + +if __name__ == "__main__": + is_ok = True + + # 1. add_field before volume is in simulation + print("Test: add_field before volume added to simulation") + orphan_box = gate.geometry.volumes.BoxVolume(name="orphan") + orphan_box.size = [10 * g4_cm, 10 * g4_cm, 10 * g4_cm] + field = fields.UniformMagneticField(name="B_orphan") + field.field_vector = [0, 0, 1 * g4_tesla] + ok = _expect_error( + lambda: orphan_box.add_field(field), + "add_field on orphan volume", + ) + is_ok = is_ok and ok + + # 2. Duplicate field on same volume + print("Test: duplicate field on same volume") + sim = gate.Simulation() + sim.world.size = [2 * g4_m, 2 * g4_m, 2 * g4_m] + sim.world.material = "G4_Galactic" + box = sim.add_volume("BoxVolume", "box") + box.size = [10 * g4_cm, 10 * g4_cm, 10 * g4_cm] + box.material = "G4_Galactic" + + f1 = fields.UniformMagneticField(name="B1") + f1.field_vector = [0, 0, 1 * g4_tesla] + box.add_field(f1) + + f2 = fields.UniformMagneticField(name="B2") + f2.field_vector = [0, 1 * g4_tesla, 0] + ok = _expect_error( + lambda: box.add_field(f2), + "second field on same volume", + ) + is_ok = is_ok and ok + + # 3. Two different fields with same name + print("Test: name collision between different field objects") + sim = gate.Simulation() + sim.world.size = [2 * g4_m, 2 * g4_m, 2 * g4_m] + sim.world.material = "G4_Galactic" + + box_a = sim.add_volume("BoxVolume", "box_a") + box_a.size = [10 * g4_cm, 10 * g4_cm, 10 * g4_cm] + box_a.material = "G4_Galactic" + box_a.translation = [-20 * g4_cm, 0, 0] + + box_b = sim.add_volume("BoxVolume", "box_b") + box_b.size = [10 * g4_cm, 10 * g4_cm, 10 * g4_cm] + box_b.material = "G4_Galactic" + box_b.translation = [20 * g4_cm, 0, 0] + + fa = fields.UniformMagneticField(name="B_same_name") + fa.field_vector = [0, 0, 1 * g4_tesla] + box_a.add_field(fa) + + fb = fields.UniformMagneticField(name="B_same_name") + fb.field_vector = [0, 1 * g4_tesla, 0] + ok = _expect_error( + lambda: box_b.add_field(fb), + "name collision", + ) + is_ok = is_ok and ok + + # 4. Same field on two volumes (should be ok) + print("Test: same field on two volumes") + sim = gate.Simulation() + sim.world.size = [2 * g4_m, 2 * g4_m, 2 * g4_m] + sim.world.material = "G4_Galactic" + + box1 = sim.add_volume("BoxVolume", "box1") + box1.size = [10 * g4_cm, 10 * g4_cm, 10 * g4_cm] + box1.material = "G4_Galactic" + box1.translation = [-20 * g4_cm, 0, 0] + + box2 = sim.add_volume("BoxVolume", "box2") + box2.size = [10 * g4_cm, 10 * g4_cm, 10 * g4_cm] + box2.material = "G4_Galactic" + box2.translation = [20 * g4_cm, 0, 0] + + shared = fields.UniformMagneticField(name="B_shared") + shared.field_vector = [0, 0, 1 * g4_tesla] + try: + box1.add_field(shared) + box2.add_field(shared) + ok = ( + shared.attached_to == ["box1", "box2"] + and box1.field == "B_shared" + and box2.field == "B_shared" + ) + print(f" {'OK' if ok else 'FAIL'} — same field on two volumes") + except Exception as e: + print(f" FAIL — unexpected error: {e}") + ok = False + is_ok = is_ok and ok + + # 5. Custom field with non-callable + print("Test: custom field with non-callable field_function") + 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(), + "non-callable field_function", + ) + is_ok = is_ok and ok + + # 6. Custom field with wrong component count + print("Test: custom B field returning 6 components instead of 3") + 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(), + "wrong component count", + ) + is_ok = is_ok and ok + + utility.test_ok(is_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 new file mode 100644 index 000000000..d2db37b87 --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_custom_vs_native_B.py @@ -0,0 +1,158 @@ +#!/usr/bin/env python3 +""" +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 +same constant vector. A proton source fires into each box. + +Checks that: + - Exit positions and energies agree within numerical tolerance. + + - Custom result agrees with the analytical cyclotron radius. +""" + +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, +) + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, output_folder="test099_fields") + + By = 5 * g4_tesla + T = 10 * g4_MeV + + sim = gate.Simulation() + sim.g4_verbose = False + sim.visu = False + sim.number_of_threads = 1 + sim.random_seed = 42 + sim.output_dir = paths.output + + world = sim.world + world.size = [3 * g4_m, 1 * g4_m, 1 * g4_m] + world.material = "G4_Galactic" + + # Box with native field + box_native = sim.add_volume("BoxVolume", "box_native") + box_native.size = [50 * g4_cm, 50 * g4_cm, 50 * g4_cm] + box_native.material = "G4_Galactic" + box_native.translation = [-80 * g4_cm, 0, 0] + + native_field = fields.UniformMagneticField(name="B_native") + native_field.field_vector = [0, By, 0] + box_native.add_field(native_field) + + # Box with custom trampoline field + box_custom = sim.add_volume("BoxVolume", "box_custom") + box_custom.size = [50 * g4_cm, 50 * g4_cm, 50 * g4_cm] + box_custom.material = "G4_Galactic" + box_custom.translation = [80 * g4_cm, 0, 0] + + def custom_uniform_B(x, y, z, t): + return [0, By, 0] + + custom_field = fields.CustomMagneticField( + name="B_custom", field_function=custom_uniform_B + ) + box_custom.add_field(custom_field) + + # Source for native box + src_native = sim.add_source("GenericSource", "src_native") + src_native.particle = "proton" + src_native.n = 1 + src_native.energy.type = "mono" + src_native.energy.mono = T + src_native.position.type = "point" + src_native.position.translation = [-80 * g4_cm, 0, -100 * g4_cm] + src_native.direction.type = "momentum" + src_native.direction.momentum = [0, 0, 1] + + # Source for custom box + src_custom = sim.add_source("GenericSource", "src_custom") + src_custom.particle = "proton" + src_custom.n = 1 + src_custom.energy.type = "mono" + src_custom.energy.mono = T + src_custom.position.type = "point" + src_custom.position.translation = [80 * g4_cm, 0, -100 * g4_cm] + src_custom.direction.type = "momentum" + src_custom.direction.momentum = [0, 0, 1] + + # Phase space actor for native box + phsp_native = sim.add_actor("PhaseSpaceActor", "phsp_native") + phsp_native.attached_to = "box_native" + phsp_native.attributes = ["PostKineticEnergy", "PostPosition"] + phsp_native.output_filename = paths.output / "test099_native_B.root" + phsp_native.steps_to_store = "exiting" + + # Phase space actor for custom box + phsp_custom = sim.add_actor("PhaseSpaceActor", "phsp_custom") + phsp_custom.attached_to = "box_custom" + phsp_custom.attributes = ["PostKineticEnergy", "PostPosition"] + phsp_custom.output_filename = paths.output / "test099_custom_B.root" + phsp_custom.steps_to_store = "exiting" + + sim.run() + + # Read results + df_native = uproot.open(str(paths.output / "test099_native_B.root"))[ + "phsp_native;1" + ].arrays(library="pd") + df_custom = uproot.open(str(paths.output / "test099_custom_B.root"))[ + "phsp_custom;1" + ].arrays(library="pd") + + # Transform to common reference frame (box centres) and compare + native_x = df_native["PostPosition_X"].values - (-80 * g4_cm) + custom_x = df_custom["PostPosition_X"].values - (80 * g4_cm) + native_y = df_native["PostPosition_Y"].values + custom_y = df_custom["PostPosition_Y"].values + native_z = df_native["PostPosition_Z"].values + custom_z = df_custom["PostPosition_Z"].values + native_KE = df_native["PostKineticEnergy"].values + custom_KE = df_custom["PostKineticEnergy"].values + + r_TOL = 1e-3 * g4_mm + e_TOL = 1e-3 * g4_MeV + + # Compare native vs custom + dx = np.abs(native_x - custom_x) + dy = np.abs(native_y - custom_y) + dz = np.abs(native_z - custom_z) + dKE = np.abs(native_KE - custom_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"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) + box_half_z = 250 * g4_mm + 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}") + + 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 new file mode 100644 index 000000000..0492b049c --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_custom_vs_native_E.py @@ -0,0 +1,156 @@ +#!/usr/bin/env python3 +""" +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 +same constant vector. A proton source fires into each box. +Exit positions and energies must agree within numerical tolerance. + +Also cross-checks the custom result against the analytical expectation +dKE = q * E * dx. +""" + +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_eplus, + PROTON_MASS, +) + +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 + + sim = gate.Simulation() + sim.g4_verbose = False + sim.visu = False + sim.number_of_threads = 1 + sim.random_seed = 42 + sim.output_dir = paths.output + + world = sim.world + world.size = [3 * g4_m, 1 * g4_m, 1 * g4_m] + world.material = "G4_Galactic" + + # Box with native field + box_native = sim.add_volume("BoxVolume", "box_native") + box_native.size = [50 * g4_cm, 50 * g4_cm, 50 * g4_cm] + box_native.material = "G4_Galactic" + box_native.translation = [-80 * g4_cm, 0, 0] + + native_field = fields.UniformElectricField(name="E_native") + native_field.field_vector = [Ex, 0, 0] + box_native.add_field(native_field) + + # Box with custom trampoline field + box_custom = sim.add_volume("BoxVolume", "box_custom") + box_custom.size = [50 * g4_cm, 50 * g4_cm, 50 * g4_cm] + box_custom.material = "G4_Galactic" + box_custom.translation = [80 * g4_cm, 0, 0] + + def custom_uniform_E(x, y, z, t): + return [Ex, 0, 0] + + custom_field = fields.CustomElectricField( + name="E_custom", field_function=custom_uniform_E + ) + box_custom.add_field(custom_field) + + # Source for native box + src_native = sim.add_source("GenericSource", "src_native") + src_native.particle = "proton" + src_native.n = 1 + src_native.energy.type = "mono" + src_native.energy.mono = T + src_native.position.type = "point" + src_native.position.translation = [-80 * g4_cm, 0, -100 * g4_cm] + src_native.direction.type = "momentum" + src_native.direction.momentum = [0, 0, 1] + + # Source for custom box + src_custom = sim.add_source("GenericSource", "src_custom") + src_custom.particle = "proton" + src_custom.n = 1 + src_custom.energy.type = "mono" + src_custom.energy.mono = T + src_custom.position.type = "point" + src_custom.position.translation = [80 * g4_cm, 0, -100 * g4_cm] + src_custom.direction.type = "momentum" + src_custom.direction.momentum = [0, 0, 1] + + # Phase space actor for native box + phsp_native = sim.add_actor("PhaseSpaceActor", "phsp_native") + phsp_native.attached_to = "box_native" + phsp_native.attributes = ["PostKineticEnergy", "PostPosition"] + phsp_native.output_filename = paths.output / "test099_native_E.root" + phsp_native.steps_to_store = "exiting" + + # Phase space actor for custom box + phsp_custom = sim.add_actor("PhaseSpaceActor", "phsp_custom") + phsp_custom.attached_to = "box_custom" + phsp_custom.attributes = ["PostKineticEnergy", "PostPosition"] + phsp_custom.output_filename = paths.output / "test099_custom_E.root" + phsp_custom.steps_to_store = "exiting" + + sim.run() + + # Read results + df_native = uproot.open(str(paths.output / "test099_native_E.root"))[ + "phsp_native;1" + ].arrays(library="pd") + df_custom = uproot.open(str(paths.output / "test099_custom_E.root"))[ + "phsp_custom;1" + ].arrays(library="pd") + + # Transform to common reference frame (box centres) and compare + native_x = df_native["PostPosition_X"].values - (-80 * g4_cm) + custom_x = df_custom["PostPosition_X"].values - (80 * g4_cm) + native_y = df_native["PostPosition_Y"].values + custom_y = df_custom["PostPosition_Y"].values + native_z = df_native["PostPosition_Z"].values + custom_z = df_custom["PostPosition_Z"].values + native_KE = df_native["PostKineticEnergy"].values + custom_KE = df_custom["PostKineticEnergy"].values + + r_TOL = 0.01 * g4_mm + e_TOL = 0.01 * g4_MeV + + # Compare native vs custom + dx = np.abs(native_x - custom_x) + dy = np.abs(native_y - custom_y) + dz = np.abs(native_z - custom_z) + dKE = np.abs(native_KE - custom_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"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}" + ) + + 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_helpers.py b/opengate/tests/src/geometry/test099_fields_helpers.py new file mode 100644 index 000000000..f47617854 --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_helpers.py @@ -0,0 +1,84 @@ +#!/usr/bin/env python3 +""" +Shared helpers for the test099_fields_* tests. +""" + +import numpy as np + +import opengate as gate +from opengate.geometry import fields + +# Units +g4_m = gate.g4_units.m +g4_cm = gate.g4_units.cm +g4_mm = gate.g4_units.mm +g4_tesla = gate.g4_units.tesla +g4_MeV = gate.g4_units.MeV +g4_volt = gate.g4_units.volt +g4_eplus = gate.g4_units.eplus + +# Proton mass +PROTON_MASS = 938.27208943 * g4_MeV + + +def cyclotron_radius(T, B, m, q): + """ + Relativistic cyclotron radius: r = p / (qBc). + """ + E_tot = T + m + p = np.sqrt(E_tot**2 - m**2) + return p / (q * B * 0.299792458) / g4_m + + +def build_field_simulation( + field_obj, + kinetic_energy=10 * gate.g4_units.MeV, + particle="proton", + n_particles=1, + phsp_output_filename="phsp.root", + n_threads=1, + seed=42, + output_dir=".", +): + """ + Build a standard simulation with: + - vacuum world (1 m^3) + - vacuum box (50 cm^3) with the given field attached + - proton source at z = -1 m, shooting along +z + - PhaseSpaceActor recording exiting particles + Returns (simulation, phsp_actor). + """ + sim = gate.Simulation() + sim.g4_verbose = False + sim.visu = False + sim.number_of_threads = n_threads + sim.random_seed = seed + sim.output_dir = output_dir + + world = sim.world + world.size = [1 * g4_m, 1 * g4_m, 1 * g4_m] + world.material = "G4_Galactic" + + box = sim.add_volume("BoxVolume", "field_box") + box.size = [50 * g4_cm, 50 * g4_cm, 50 * g4_cm] + box.material = "G4_Galactic" + box.add_field(field_obj) + + source = sim.add_source("GenericSource", "particle_source") + source.particle = particle + source.n = n_particles + source.energy.type = "mono" + source.energy.mono = kinetic_energy + source.position.type = "point" + source.position.translation = [0, 0, -100 * g4_cm] + source.direction.type = "momentum" + source.direction.momentum = [0, 0, 1] + + phsp = sim.add_actor("PhaseSpaceActor", "phsp") + phsp.attached_to = box.name + phsp.attributes = ["PostKineticEnergy", "PostPosition"] + phsp.output_filename = phsp_output_filename + phsp.steps_to_store = "exiting" + phsp.root_output.write_to_disk = True + + return sim, phsp diff --git a/opengate/tests/src/geometry/test099_fields_serialization.py b/opengate/tests/src/geometry/test099_fields_serialization.py new file mode 100644 index 000000000..727bfc09f --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_serialization.py @@ -0,0 +1,151 @@ +#!/usr/bin/env python3 +""" +Test 099 — Field serialization round-trip. + +Tests that non-custom field types survive serialization +and that custom fields correctly refuse it. +""" + +import opengate as gate +from opengate.geometry import fields +from opengate.tests import utility + +from test099_fields_helpers import g4_tesla, g4_m, g4_cm, g4_mm, g4_volt + + +def _make_sim_with_box(): + """Create a minimal simulation with one box volume.""" + sim = gate.Simulation() + sim.number_of_threads = 1 + world = sim.world + world.size = [2 * g4_m, 2 * g4_m, 2 * g4_m] + world.material = "G4_Galactic" + + box = sim.add_volume("BoxVolume", "box") + box.size = [50 * g4_cm, 50 * g4_cm, 50 * g4_cm] + box.material = "G4_Galactic" + return sim, box + + +if __name__ == "__main__": + is_ok = True + + # Uniform magnetic field round-trip + sim, box = _make_sim_with_box() + field = fields.UniformMagneticField(name="B_uniform") + field.field_vector = [0, 1 * g4_tesla, 0] + field.delta_chord = 0.5 * g4_mm + box.add_field(field) + + d = sim.to_dictionary() + assert "fields" in d["volume_manager"] + assert "B_uniform" in d["volume_manager"]["fields"] + fd = d["volume_manager"]["fields"]["B_uniform"] + assert fd["object_type"] == "UniformMagneticField" + assert fd["attached_to"] == ["box"] + assert fd["user_info"]["field_vector"] == [0, 1 * g4_tesla, 0] + + sim2 = gate.Simulation() + sim2.from_dictionary(d) + restored = sim2.volume_manager.fields["B_uniform"] + ok = ( + isinstance(restored, fields.UniformMagneticField) + and restored.field_vector == [0, 1 * g4_tesla, 0] + and restored.delta_chord == 0.5 * g4_mm + and restored.attached_to == ["box"] + and sim2.volume_manager.volumes["box"].field == "B_uniform" + ) + is_ok = is_ok and ok + print(f"Uniform B round-trip: {'OK' if ok else 'FAIL'}") + + # Quadrupole round-trip + sim, box = _make_sim_with_box() + field = fields.QuadrupoleMagneticField(name="B_quad") + field.gradient = 10 * 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_quad"] + ok = ( + isinstance(restored, fields.QuadrupoleMagneticField) + and restored.gradient == 10 * g4_tesla / g4_m + and restored.attached_to == ["box"] + ) + is_ok = is_ok and ok + print(f"Quadrupole round-trip: {'OK' if ok else 'FAIL'}") + + # Uniform electric field round-trip + sim, box = _make_sim_with_box() + field = fields.UniformElectricField(name="E_uniform") + field.field_vector = [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_uniform"] + ok = ( + isinstance(restored, fields.UniformElectricField) + and restored.field_vector == [0, 0, 1e6 * g4_volt / g4_m] + and restored.attached_to == ["box"] + ) + is_ok = is_ok and ok + print(f"Uniform E round-trip: {'OK' if ok else 'FAIL'}") + + # Custom fields refuse serialization + field = fields.CustomMagneticField(name="B_custom") + field.field_function = lambda x, y, z, t: [0, 0, 0] + try: + field.to_dictionary() + ok = False # should not reach here + except NotImplementedError: + ok = True + is_ok = is_ok and ok + print(f"Custom B serialization raises: {'OK' if ok else 'FAIL'}") + + field = fields.CustomElectricField(name="E_custom") + field.field_function = lambda x, y, z, t: [0, 0, 0] + try: + field.to_dictionary() + ok = False + except NotImplementedError: + ok = True + is_ok = is_ok and ok + print(f"Custom E serialization raises: {'OK' if ok else 'FAIL'}") + + # Multiple volume attachment round-trip + sim = gate.Simulation() + sim.number_of_threads = 1 + sim.world.size = [2 * g4_m, 2 * g4_m, 2 * g4_m] + sim.world.material = "G4_Galactic" + + box1 = sim.add_volume("BoxVolume", "box1") + box1.size = [20 * g4_cm, 20 * g4_cm, 20 * g4_cm] + box1.material = "G4_Galactic" + box1.translation = [-30 * g4_cm, 0, 0] + + box2 = sim.add_volume("BoxVolume", "box2") + box2.size = [20 * g4_cm, 20 * g4_cm, 20 * g4_cm] + box2.material = "G4_Galactic" + box2.translation = [30 * g4_cm, 0, 0] + + field = fields.UniformMagneticField(name="B_shared") + field.field_vector = [0, 0, 2 * g4_tesla] + box1.add_field(field) + box2.add_field(field) + + d = sim.to_dictionary() + sim2 = gate.Simulation() + sim2.from_dictionary(d) + restored = sim2.volume_manager.fields["B_shared"] + ok = ( + set(restored.attached_to) == {"box1", "box2"} + and sim2.volume_manager.volumes["box1"].field == "B_shared" + and sim2.volume_manager.volumes["box2"].field == "B_shared" + ) + is_ok = is_ok and ok + print(f"Multi-volume round-trip: {'OK' if ok else 'FAIL'}") + + utility.test_ok(is_ok)