EM Field Enhancement: local coordinates, dynamic placements and mapped fields

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 <pybind11/pybind11.h>
+
+namespace py = pybind11;
+
+#include "G4MagneticField.hh"
+#include "G4SextupoleMagField.hh"
+
+void init_G4SextupoleMagField(py::module &m) {
+
+  py::class_<G4SextupoleMagField, G4MagneticField,
+             std::unique_ptr<G4SextupoleMagField, py::nodelete>>(
+      m, "G4SextupoleMagField")
+
+      .def(py::init<G4double>())
+
+      ;
+}

core/opengate_core/opengate_core.cpp

@@ -187,6 +187,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 &);
@@ -209,6 +211,14 @@ 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 &);
 
@@ -549,6 +559,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);
@@ -560,6 +571,13 @@ PYBIND11_MODULE(opengate_core, 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);
@@ -668,6 +686,7 @@ PYBIND11_MODULE(opengate_core, m) {
   init_GateRunAction(m);
   init_GateEventAction(m);
   init_GateTrackingAction(m);
+  // init_GateField(m);
 
   init_GateDoseActor(m);
   init_GateTLEDoseActor(m);

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<G4ThreeVector> translations,
+    std::vector<G4RotationMatrix> 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);
+}

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 <vector>
+
+class G4VSolid;
+
+// GATE wrapper for G4ElectroMagneticField.
+class GateElectroMagneticField : public G4ElectroMagneticField,
+                                 public GateFieldBase {
+public:
+  // constructor
+  GateElectroMagneticField(G4ElectroMagneticField *inner, const G4VSolid *solid,
+                           std::vector<G4ThreeVector> translations,
+                           std::vector<G4RotationMatrix> 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

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 <cmath>
+#include <limits>
+#include <sstream>
+#include <stdexcept>
+
+// constructor
+GateFieldBase::GateFieldBase(const G4VSolid *solid,
+                             std::vector<G4ThreeVector> translations,
+                             std::vector<G4RotationMatrix> 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<G4ThreeVector> translations,
+                                  std::vector<G4RotationMatrix> 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<double>::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);
+}

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 <vector>
+
+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<G4ThreeVector> translations,
+                std::vector<G4RotationMatrix> rotations, double deltaChordMM);
+
+  // update the transforms (e.g. after a geometry change between runs)
+  void SetTransforms(std::vector<G4ThreeVector> translations,
+                     std::vector<G4RotationMatrix> 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 <typename LocalFieldFn>
+
+  // 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<G4AffineTransform> m_transforms;
+
+  static constexpr double kMaxCurvatureRadiusMM = 100.0e3;
+};
+
+#endif // GateFieldBase_h