Skip to content

Implement electric and magnetic fields in GATE #1003

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
tbaudier wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Implement electric and magnetic fields in GATE #1003

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 26 additions & 0 deletions core/opengate_core/g4_bindings/pyG4ChordFinder.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 <pybind11/pybind11.h>

namespace py = pybind11;

#include "G4ChordFinder.hh"
#include "G4MagIntegratorStepper.hh"
#include "G4MagneticField.hh"
#include "G4VIntegrationDriver.hh"

void init_G4ChordFinder(py::module &m) {
py::class_<G4ChordFinder, std::unique_ptr<G4ChordFinder, py::nodelete>>(
m, "G4ChordFinder")

.def(py::init<G4VIntegrationDriver *>())
.def(py::init<G4MagneticField *, G4double, G4MagIntegratorStepper *,
G4int>())

.def("SetDeltaChord", &G4ChordFinder::SetDeltaChord);
}
24 changes: 24 additions & 0 deletions core/opengate_core/g4_bindings/pyG4ClassicalRK4.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 "G4ClassicalRK4.hh"
#include "G4EquationOfMotion.hh"
#include "G4MagErrorStepper.hh"

void init_G4ClassicalRK4(py::module &m) {
// G4ClassicalRK4 inherits from G4MagErrorStepper
py::class_<G4ClassicalRK4, G4MagErrorStepper,
std::unique_ptr<G4ClassicalRK4, py::nodelete>>(m, "G4ClassicalRK4")

.def(py::init<G4EquationOfMotion *, G4int>())

;
}
75 changes: 75 additions & 0 deletions core/opengate_core/g4_bindings/pyG4ElectricField.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,75 @@
/* --------------------------------------------------
Copyright (C): OpenGATE Collaboration
This software is distributed under the terms
of the GNU Lesser General Public Licence (LGPL)
See LICENSE.md for further details
-------------------------------------------------- */

#include <pybind11/pybind11.h>
#include <pybind11/stl.h>

namespace py = pybind11;

#include "G4ElectricField.hh"
#include "G4ElectroMagneticField.hh"

// Trampoline class to allow Python to override GetFieldValue for electric
// fields. Inherits from G4ElectricField which already implements
// DoesFieldChangeEnergy() = true.
class PyG4ElectricField : public G4ElectricField {
public:
using G4ElectricField::G4ElectricField;

void GetFieldValue(const G4double Point[4], G4double *field) const override {
// Always initialize output to a safe value: [Bx, By, Bz, Ex, Ey, Ez].
field[0] = 0.;
field[1] = 0.;
field[2] = 0.;
field[3] = 0.;
field[4] = 0.;
field[5] = 0.;

py::gil_scoped_acquire gil;

// Convert Point to Python list
py::list pyPoint;
pyPoint.append(Point[0]);
pyPoint.append(Point[1]);
pyPoint.append(Point[2]);
pyPoint.append(Point[3]);

// Get the Python override
py::function override = py::get_override(this, "GetFieldValue");
if (override) {
py::object result = override(pyPoint);

if (!result.is_none()) {
py::sequence field_seq = result.cast<py::sequence>();

if (py::len(field_seq) != 3) {
throw std::invalid_argument(
"GetFieldValue for G4ElectricField must return exactly 3 "
"components [Ex, Ey, Ez]");
}

// User returned [Ex, Ey, Ez]
field[3] = field_seq[0].cast<G4double>();
field[4] = field_seq[1].cast<G4double>();
field[5] = field_seq[2].cast<G4double>();
}
}
}
};

void init_G4ElectricField(py::module &m) {

py::class_<G4ElectricField, G4ElectroMagneticField, PyG4ElectricField,
std::unique_ptr<G4ElectricField, py::nodelete>>(m,
"G4ElectricField")

.def(py::init<>())

.def("DoesFieldChangeEnergy", &G4ElectricField::DoesFieldChangeEnergy)

;
}
80 changes: 80 additions & 0 deletions core/opengate_core/g4_bindings/pyG4ElectroMagneticField.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@
/* --------------------------------------------------
Copyright (C): OpenGATE Collaboration
This software is distributed under the terms
of the GNU Lesser General Public Licence (LGPL)
See LICENSE.md for further details
-------------------------------------------------- */

#include <pybind11/pybind11.h>
#include <pybind11/stl.h>

namespace py = pybind11;

#include "G4ElectroMagneticField.hh"
#include "G4Field.hh"

// Trampoline class to allow Python to override GetFieldValue for
// electromagnetic fields.
class PyG4ElectroMagneticField : public G4ElectroMagneticField {
public:
using G4ElectroMagneticField::G4ElectroMagneticField;

void GetFieldValue(const G4double Point[4], G4double *field) const override {
// Always initialize output to a safe value: [Bx, By, Bz, Ex, Ey, Ez].
field[0] = 0.;
field[1] = 0.;
field[2] = 0.;
field[3] = 0.;
field[4] = 0.;
field[5] = 0.;

py::gil_scoped_acquire gil;

// Convert Point to Python list
py::list pyPoint;
pyPoint.append(Point[0]);
pyPoint.append(Point[1]);
pyPoint.append(Point[2]);
pyPoint.append(Point[3]);

// Try to get the Python override
py::function override = py::get_override(this, "GetFieldValue");
if (override) {
// Call Python implementation and expect [Bx, By, Bz, Ex, Ey, Ez]
py::object result = override(pyPoint);
if (!result.is_none()) {
py::sequence field_seq = result.cast<py::sequence>();
size_t n = py::len(field_seq);

if (n != 6) {
throw std::invalid_argument(
"GetFieldValue for G4ElectroMagneticField must return exactly 6 "
"components [Bx, By, Bz, Ex, Ey, Ez]");
}

for (size_t i = 0; i < n && i < 6; ++i) {
field[i] = field_seq[i].cast<G4double>();
}
}
}
}

G4bool DoesFieldChangeEnergy() const override {
PYBIND11_OVERRIDE_PURE(G4bool, G4ElectroMagneticField,
DoesFieldChangeEnergy);
}
};

void init_G4ElectroMagneticField(py::module &m) {

py::class_<G4ElectroMagneticField, G4Field, PyG4ElectroMagneticField,
std::unique_ptr<G4ElectroMagneticField, py::nodelete>>(
m, "G4ElectroMagneticField")

.def(py::init<>())

.def("DoesFieldChangeEnergy",
&G4ElectroMagneticField::DoesFieldChangeEnergy)

;
}
28 changes: 28 additions & 0 deletions core/opengate_core/g4_bindings/pyG4EqMagElectricField.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 <pybind11/pybind11.h>

namespace py = pybind11;

#include "G4ElectroMagneticField.hh"
#include "G4EqMagElectricField.hh"
#include "G4EquationOfMotion.hh"

void init_G4EqMagElectricField(py::module &m) {

py::class_<G4EqMagElectricField, G4EquationOfMotion,
std::unique_ptr<G4EqMagElectricField, py::nodelete>>(
m, "G4EqMagElectricField")

.def(py::init<G4ElectroMagneticField *>())

.def("SetChargeMomentumMass",
&G4EqMagElectricField::SetChargeMomentumMass)

;
}
23 changes: 23 additions & 0 deletions core/opengate_core/g4_bindings/pyG4EquationOfMotion.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
/* --------------------------------------------------
Copyright (C): OpenGATE Collaboration
This software is distributed under the terms
of the GNU Lesser General Public Licence (LGPL)
See LICENSE.md for further details
-------------------------------------------------- */

#include <pybind11/pybind11.h>

namespace py = pybind11;

#include "G4EquationOfMotion.hh"

void init_G4EquationOfMotion(py::module &m) {

py::class_<G4EquationOfMotion,
std::unique_ptr<G4EquationOfMotion, py::nodelete>>(
m, "G4EquationOfMotion")

.def("GetFieldObj",
py::overload_cast<>(&G4EquationOfMotion::GetFieldObj, py::const_),
py::return_value_policy::reference_internal);
}
80 changes: 80 additions & 0 deletions core/opengate_core/g4_bindings/pyG4Field.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@
/* --------------------------------------------------
Copyright (C): OpenGATE Collaboration
This software is distributed under the terms
of the GNU Lesser General Public Licence (LGPL)
See LICENSE.md for further details
-------------------------------------------------- */

#include <pybind11/pybind11.h>
#include <pybind11/stl.h>

namespace py = pybind11;

#include "G4Field.hh"

/*
* Trampoline class to allow Python to override GetFieldValue and
* DoesFieldChangeEnergy. This enables users to define custom fields in Python.
*
* The GetFieldValue method receives a point (x, y, z, t) and must return
* the field components at that point. The number of components depends on
* the field type:
* - Magnetic field: 3 components (Bx, By, Bz)
* - Electric field: 3 components (Ex, Ey, Ez) at indices 3-5
* - Electromagnetic field: 6 components (Bx, By, Bz, Ex, Ey, Ez)
*/
class PyG4Field : public G4Field {
public:
// Inherit the constructors
using G4Field::G4Field;

void GetFieldValue(const G4double Point[4],
G4double *fieldArr) const override {
py::gil_scoped_acquire gil;
// Always initialize output to a safe value.
for (size_t i = 0; i < G4Field::MAX_NUMBER_OF_COMPONENTS; ++i) {
fieldArr[i] = 0.;
}

// Convert Point to Python list
py::list pyPoint;
pyPoint.append(Point[0]);
pyPoint.append(Point[1]);
pyPoint.append(Point[2]);
pyPoint.append(Point[3]);

// Try to get the Python override
py::function override = py::get_override(this, "GetFieldValue");
if (override) {
// Call Python implementation and expect a list back
py::object result = override(pyPoint);
if (!result.is_none()) {
py::list field_list = result.cast<py::list>();
size_t n = py::len(field_list);
for (size_t i = 0; i < n && i < G4Field::MAX_NUMBER_OF_COMPONENTS;
++i) {
fieldArr[i] = field_list[i].cast<G4double>();
}
}
}
}

G4bool DoesFieldChangeEnergy() const override {
PYBIND11_OVERRIDE_PURE(G4bool, G4Field, DoesFieldChangeEnergy);
}
};

void init_G4Field(py::module &m) {

py::class_<G4Field, PyG4Field, std::unique_ptr<G4Field, py::nodelete>>(
m, "G4Field")

.def(py::init<G4bool>())
.def(py::init<>())

.def("DoesFieldChangeEnergy", &G4Field::DoesFieldChangeEnergy)
.def("IsGravityActive", &G4Field::IsGravityActive)
.def("SetGravityActive", &G4Field::SetGravityActive)

;
}
Loading
Loading