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/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 8eb2cfa39b..f11b571323 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 &); @@ -204,12 +206,36 @@ 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 &); void init_G4ChordFinder(py::module &); +void init_GateMagneticField(py::module &); +void init_GateElectroMagneticField(py::module &); +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 &); + // geometry/solids void init_G4Box(py::module &); @@ -573,6 +599,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); @@ -581,9 +608,24 @@ 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); + init_GateMagneticField(m); + init_GateElectroMagneticField(m); + init_GateUniformElectroMagneticField(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/GateElectroMagneticField.cpp b/core/opengate_core/opengate_lib/GateElectroMagneticField.cpp new file mode 100644 index 0000000000..fc143ca6ad --- /dev/null +++ b/core/opengate_core/opengate_lib/GateElectroMagneticField.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 "GateElectroMagneticField.h" + +GateElectroMagneticField::GateElectroMagneticField( + G4ElectroMagneticField *inner, const G4VSolid *solid, + std::vector translations, + std::vector rotations, double deltaChordMM) + : G4ElectroMagneticField(), + GateFieldBase(solid, std::move(translations), std::move(rotations), + deltaChordMM), + m_inner(inner) {} + +void GateElectroMagneticField::GetFieldValue(const G4double Point[4], + G4double *BEfield) 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); + }; + + applyTransforms(Point, BEfield, 6, localFieldFunc); +} diff --git a/core/opengate_core/opengate_lib/GateElectroMagneticField.h b/core/opengate_core/opengate_lib/GateElectroMagneticField.h new file mode 100644 index 0000000000..10a3e22c0a --- /dev/null +++ b/core/opengate_core/opengate_lib/GateElectroMagneticField.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 GateElectroMagneticField_h +#define GateElectroMagneticField_h + +#include "G4ElectroMagneticField.hh" +#include "GateFieldBase.h" +#include + +class G4VSolid; + +// GATE wrapper for G4ElectroMagneticField. +class GateElectroMagneticField : public G4ElectroMagneticField, + 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; } + +private: + // the inner field in the volume's local frame + G4ElectroMagneticField *m_inner; +}; + +#endif // GateElectroMagneticField_h diff --git a/core/opengate_core/opengate_lib/GateFieldBase.cpp b/core/opengate_core/opengate_lib/GateFieldBase.cpp new file mode 100644 index 0000000000..f4893fe0d4 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateFieldBase.cpp @@ -0,0 +1,111 @@ +/* -------------------------------------------------- + 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 "GateFieldBase.h" + +#include "G4VSolid.hh" +#include "globals.hh" + +#include +#include +#include +#include + +// constructor +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)) { + if (solid == nullptr) + throw std::invalid_argument("GateFieldBase: solid must not be null"); + + if (translations.size() != rotations.size() || translations.empty()) + throw std::invalid_argument("GateFieldBase: translations and rotations " + "must be non-empty and have the same size"); + + 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 GateFieldBase::SetTransforms(std::vector translations, + std::vector rotations) { + if (translations.size() != rotations.size() || translations.empty()) + throw std::invalid_argument( + "GateFieldBase::SetTransforms: translations and rotations must be " + "non-empty and have the same size"); + + 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 +G4ThreeVector GateFieldBase::findContainingPlacement( + const G4ThreeVector &worldPoint, + const G4AffineTransform *&outTransform) const { + 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; + } + + const double d = m_solid->DistanceToIn(localPoint); + if (d < minDistToSurface) { + minDistToSurface = d; + closestIdx = i; + closestLocal = localPoint; + } + } + + // 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 << "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: " + << m_fallbackFatalDistanceMM << " mm.\n" + << " This likely indicates a real bug in the geometry or field " + "setup.\n"; + + 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 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/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/GateMagneticField.cpp b/core/opengate_core/opengate_lib/GateMagneticField.cpp new file mode 100644 index 0000000000..d2e38e71e9 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateMagneticField.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 "GateMagneticField.h" + +// constructor +GateMagneticField::GateMagneticField(G4MagneticField *inner, + const G4VSolid *solid, + std::vector translations, + std::vector rotations, + double deltaChordMM) + : G4MagneticField(), GateFieldBase(solid, std::move(translations), + std::move(rotations), deltaChordMM), + m_inner(inner) {} + +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); + }; + + applyTransforms(Point, Bfield, 3, localFieldFunc); +} diff --git a/core/opengate_core/opengate_lib/GateMagneticField.h b/core/opengate_core/opengate_lib/GateMagneticField.h new file mode 100644 index 0000000000..e05b3432b1 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateMagneticField.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 GateMagneticField_h +#define GateMagneticField_h + +#include "G4MagneticField.hh" +#include "GateFieldBase.h" +#include + +class G4VSolid; + +// GATE wrapper for G4MagneticField +class GateMagneticField : public G4MagneticField, public GateFieldBase { +public: + // constructor + GateMagneticField(G4MagneticField *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 *Bfield) const override; + +private: + // the inner field in the volume's local frame + G4MagneticField *m_inner; +}; + +#endif // GateMagneticField_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/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/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/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/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/pyGateMagneticField.cpp b/core/opengate_core/opengate_lib/pyGateMagneticField.cpp new file mode 100644 index 0000000000..c1a26b588c --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateMagneticField.cpp @@ -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 + -------------------------------------------------- */ + +#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."); +} 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/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/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/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/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 d4de9a5a7e..e56e75bf77 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 @@ -64,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 ~~~~~~~~~~~~~~~ @@ -76,26 +112,50 @@ 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 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 ~~~~~~~~~~~~~~~~~~~~~~ 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_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]``. .. code-block:: python @@ -107,8 +167,21 @@ 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 +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. @@ -120,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. @@ -145,15 +221,40 @@ 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 ---------------------------- -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 +278,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 +309,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,11 +318,11 @@ 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. -**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 @@ -243,8 +349,18 @@ 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_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. +- ``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 @@ -260,16 +376,30 @@ Class reference .. autoclass:: opengate.geometry.fields.QuadrupoleMagneticField :members: +.. autoclass:: opengate.geometry.fields.SextupoleMagneticField + :members: + .. 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: .. automethod:: opengate.geometry.volumes.VolumeBase.add_field 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..a7aff20502 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -953,7 +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) - volume_obj.g4_field_manager = field.create_field_manager() + 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 a487bcb91f..d490fca4a2 100644 --- a/opengate/geometry/fields.py +++ b/opengate/geometry/fields.py @@ -1,14 +1,19 @@ +from typing import Any + +import numpy as np 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) -# ! - 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. @@ -19,6 +24,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 @@ -27,6 +33,18 @@ class FieldBase(GateObject): max_epsilon_step: float user_info_defaults = { + "stepper": ( + "DormandPrince745", + { + "doc": ( + "Integration stepper type. " + "General-purpose (any field): 'DormandPrince745' (default), 'ClassicalRK4', 'CashKarpRKF45', " + "'BogackiShampine45', 'BogackiShampine23', 'DormandPrinceRK56', " + "'DormandPrinceRK78'. " + "Magnetic-only: 'NystromRK4', 'ExactHelixStepper'." + ), + }, + ), "step_minimum": ( 1e-2 * g4_units.mm, { @@ -70,106 +88,176 @@ def __init__(self, *args, **kwargs) -> None: self.g4_field = None self._g4_runtime_objects = [] + self._field_volume_obj: Any = 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).""" return self._field_changes_energy - def create_field_manager(self) -> g4.G4FieldManager: + def refresh_transforms(self) -> None: + """Recompute and push cached world-to-local transforms after dynamic + geometry changes. + """ + 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()) + 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]: + # 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_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)) + g4_rotations.append(rot_np_as_g4(R)) + + g4_field.SetTransforms(g4_translations, g4_rotations) + + 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) - 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 _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 + ) + return ( + [vec_np_as_g4(t) for t in translations_np], + [rot_np_as_g4(r) for r in rotations_np], + ) - def _create_field(self) -> None: - """Create the G4 field object. Override.""" - msg = "_create_field() must be implemented in subclasses." - raise NotImplementedError(msg) + 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"Stepper '{self.stepper}' only supports pure magnetic fields. " + "Choose a general-purpose stepper for electric or EM fields." + ) - 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) - ) + 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) + 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, - 6, # number of variables for magnetic field = 6 (x,y,z + px,py,pz) - 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) - # Create and configure field manager 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) 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. self._g4_runtime_objects.append( { - "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, "chord_finder": self.g4_chord_finder, "field_manager": fm, + "volume": volume_obj, } ) 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, 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( + 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, volume_obj + ) + + 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,17 +266,14 @@ 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).", }, ), } - 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): @@ -206,12 +291,27 @@ class QuadrupoleMagneticField(MagneticField): ), } - def __init__(self, *args, **kwargs) -> None: - super().__init__(*args, **kwargs) + def _create_inner_field(self): + return g4.G4QuadrupoleMagField(self.gradient) + + +class SextupoleMagneticField(MagneticField): + """Sextupole magnetic field with gradient.""" - def _create_field(self) -> None: - """Create the quadrupole magnetic field.""" - self.g4_field = g4.G4QuadrupoleMagField(self.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): @@ -224,30 +324,20 @@ 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.", }, ), } - def __init__(self, *args, **kwargs) -> None: - super().__init__(*args, **kwargs) + def _create_inner_field(self): + """Create the custom magnetic field using the Python trampoline. - 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]" - ) + 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. + """ + _validate_field_function(self.field_function, "CustomMagneticField", 3) - # Create a custom G4MagneticField subclass that calls our Python function class _PyMagneticField(g4.G4MagneticField): def __init__(inner_self, callback): super().__init__() @@ -256,7 +346,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( @@ -265,82 +355,43 @@ 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_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, + 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." ) - 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 + 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( + inner, + self._field_volume_obj.g4_solid, + g4_translations, + g4_rotations, + self.delta_chord, ) - 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, - } + # 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, volume_obj ) - 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 + # _field_changes_energy is already True from ElectroMagneticField class UniformElectricField(ElectricField): @@ -353,17 +404,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): @@ -376,30 +424,14 @@ 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.", }, ), } - def __init__(self, *args, **kwargs) -> None: - super().__init__(*args, **kwargs) + def _create_inner_field(self): + _validate_field_function(self.field_function, "CustomElectricField", 3) - 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__() @@ -411,7 +443,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( @@ -419,60 +451,33 @@ 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.""" # 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 Geant4 units.", + "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 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.field_vector_E), + g4.G4ThreeVector(*self.field_vector_B), + ) class CustomElectroMagneticField(ElectroMagneticField): @@ -485,32 +490,14 @@ 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.", }, ), } - def __init__(self, *args, **kwargs) -> None: - super().__init__(*args, **kwargs) + def _create_inner_field(self): + _validate_field_function(self.field_function, "CustomElectroMagneticField", 6) - 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__() @@ -522,7 +509,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( @@ -530,24 +517,339 @@ 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. + + 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 + + +# 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, +} + +_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 with separate B and E grids on regular 3D Cartesian grids. + + 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: + 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 = { "UniformMagneticField": UniformMagneticField, "QuadrupoleMagneticField": QuadrupoleMagneticField, + "SextupoleMagneticField": SextupoleMagneticField, "CustomMagneticField": CustomMagneticField, + "MappedMagneticField": MappedMagneticField, "UniformElectricField": UniformElectricField, "CustomElectricField": CustomElectricField, + "MappedElectricField": MappedElectricField, "UniformElectroMagneticField": UniformElectroMagneticField, "CustomElectroMagneticField": CustomElectroMagneticField, + "MappedElectroMagneticField": MappedElectroMagneticField, } process_cls(FieldBase) process_cls(MagneticField) process_cls(UniformMagneticField) process_cls(QuadrupoleMagneticField) +process_cls(SextupoleMagneticField) process_cls(CustomMagneticField) +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) 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 294b77e1a6..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 @@ -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_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_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_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 new file mode 100644 index 0000000000..c094cd8b7f --- /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..36bc488c62 --- /dev/null +++ b/opengate/tests/src/geometry/test099_fields_mapped_vs_uniform_E.py @@ -0,0 +1,161 @@ +#!/usr/bin/env python3 +""" +Test 099 - MappedElectricField with uniform grid vs UniformElectricField. + +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: + - 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) 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..92e70b683e 100644 --- a/opengate/tests/src/geometry/test099_fields_serialization.py +++ b/opengate/tests/src/geometry/test099_fields_serialization.py @@ -1,11 +1,15 @@ #!/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. """ +import itertools + +import numpy as np + import opengate as gate from opengate.geometry import fields from opengate.tests import utility @@ -76,6 +80,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") @@ -148,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) 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..95a57749da --- /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-2 * g4_mm + e_TOL = 1e-2 * 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..1429bf94e5 --- /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 = False + +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)