Added the AMDM actor as described in Parisi et al 10.1002/acm2.14049

core/opengate_core/opengate_core.cpp

@@ -412,6 +412,7 @@ void init_GateVoxelizedPromptGammaAnalogActor(py::module &m);
 void init_GateFluenceActor(py::module &m);
 
 void init_GateLETActor(py::module &m);
+void init_GateAMDMActor(py::module &m);
 
 void init_GateProductionAndStoppingActor(py::module &m);
 
@@ -823,6 +824,7 @@ PYBIND11_MODULE(opengate_core, m) {
   init_GateVoxelizedPromptGammaAnalogActor(m);
   init_GateFluenceActor(m);
   init_GateLETActor(m);
+  init_GateAMDMActor(m);
   init_GateProductionAndStoppingActor(m);
   init_GateBeamQualityActor(m);
   init_GateEmCalculatorActor(m);

core/opengate_core/opengate_lib/GateAMDMActor.cpp

@@ -0,0 +1,140 @@
+/* --------------------------------------------------
+   Copyright (C): OpenGATE Collaboration
+   LGPL, see LICENSE.md.
+   Original AMDM implementation: Hermann Fuchs, April 2025.
+   -------------------------------------------------- */
+#include "GateAMDMActor.h"
+#include "GateHelpersDict.h"
+#include "GateHelpersGeometry.h"
+#include "GateHelpersImage.h"
+#include <G4RandomTools.hh>
+#include <cmath>
+#include <limits>
+#include <stdexcept>
+
+GateAMDMActor::GateAMDMActor(py::dict &user_info)
+    : GateVActor(user_info, true) {}
+
+void GateAMDMActor::InitializeUserInfo(py::dict &user_info) {
+  GateVActor::InitializeUserInfo(user_info);
+  fBins = DictGetInt(user_info, "AMDM_Bins");
+  fHitType = DictGetStr(user_info, "hit_type");
+  if (fHitType != "random" && fHitType != "pre" && fHitType != "post" &&
+      fHitType != "middle")
+    throw std::invalid_argument(
+        "AMDM hit_type must be random, pre, post or middle");
+  fTranslation = DictGetG4ThreeVector(user_info, "translation");
+  fRotation = DictGetG4RotationMatrix(user_info, "rotation");
+  fTable.Read(DictGetStr(user_info, "LUTfilename"), fBins);
+}
+
+void GateAMDMActor::InitializeCpp() {
+  GateVActor::InitializeCpp();
+  cpp_amdm_restricted_edep_image = ImageType::New();
+  cpp_amdm_delta_image = ImageType4D::New();
+  cpp_amdm_gamma_image = ImageType4D::New();
+}
+
+void GateAMDMActor::BeginOfRunActionMasterThread(int) {
+  NbOfEvent = 0;
+  // Python supplies fresh zero-filled raw buffers for each run. Only attach the
+  // 3D domain: the generic attachment helper initializes only three axes.
+  G4ThreeVector volumeTranslation;
+  G4RotationMatrix volumeRotation;
+  ComputeTransformationFromVolumeToWorld(fPhysicalVolumeName, volumeTranslation,
+                                         volumeRotation, true);
+  G4ThreeVector center;
+  const auto size =
+      cpp_amdm_restricted_edep_image->GetLargestPossibleRegion().GetSize();
+  const auto spacing = cpp_amdm_restricted_edep_image->GetSpacing();
+  for (unsigned int i = 0; i < 3; ++i)
+    center[i] = -(size[i] - 1.0) * spacing[i] / 2.0;
+  const auto spatialOrigin =
+      volumeTranslation + volumeRotation * (fRotation * center + fTranslation);
+  const auto spatialRotation = volumeRotation * fRotation;
+  ImageType::PointType origin3;
+  ImageType::DirectionType direction3;
+  for (unsigned int i = 0; i < 3; ++i) {
+    origin3[i] = spatialOrigin[i];
+    for (unsigned int j = 0; j < 3; ++j)
+      direction3[i][j] = spatialRotation(i, j);
+  }
+  cpp_amdm_restricted_edep_image->SetOrigin(origin3);
+  cpp_amdm_restricted_edep_image->SetDirection(direction3);
+  ImageType4D::PointType origin;
+  ImageType4D::DirectionType direction;
+  direction.SetIdentity();
+  origin.Fill(0);
+  for (unsigned int i = 0; i < 3; ++i) {
+    origin[i] = cpp_amdm_restricted_edep_image->GetOrigin()[i];
+    for (unsigned int j = 0; j < 3; ++j)
+      direction[i][j] = cpp_amdm_restricted_edep_image->GetDirection()[i][j];
+  }
+  for (auto image : {cpp_amdm_delta_image, cpp_amdm_gamma_image}) {
+    image->SetOrigin(origin);
+    image->SetDirection(direction);
+  }
+}
+
+void GateAMDMActor::BeginOfEventAction(const G4Event *) {
+  std::lock_guard<std::mutex> lock(fScoringMutex);
+  ++NbOfEvent;
+}
+
+std::vector<double> GateAMDMActor::Lookup(int charge, double energy) const {
+  std::vector<double> values;
+  fTable.Find(charge, energy, values);
+  return values; // empty means this charge/sample does not contribute
+}
+
+void GateAMDMActor::SteppingAction(G4Step *step) {
+  const auto pre = step->GetPreStepPoint()->GetPosition();
+  const auto post = step->GetPostStepPoint()->GetPosition();
+  auto position = post;
+  if (fHitType == "pre")
+    position = pre;
+  else if (fHitType == "middle")
+    position = pre + 0.5 * (post - pre);
+  else if (fHitType == "random")
+    position = pre + G4UniformRand() * (post - pre);
+  // As in DoseActor, the attached image domain is in world coordinates.
+  ImageType::PointType point;
+  for (unsigned int i = 0; i < 3; ++i)
+    point[i] = position[i];
+  ImageType::IndexType index;
+  if (!cpp_amdm_restricted_edep_image->TransformPhysicalPointToIndex(point,
+                                                                     index))
+    return;
+  const auto *track = step->GetTrack();
+  const int charge = static_cast<int>(track->GetDefinition()->GetPDGCharge());
+  if (charge <= 0)
+    return;
+  const int mass = track->GetDefinition()->GetBaryonNumber();
+  // GetKineticEnergy is the track's post-step sample at this scoring callback.
+  const double kineticEnergy = track->GetKineticEnergy() / CLHEP::MeV;
+  // Preserve legacy applicability, including its zero-baryon rule: positive
+  // energy / 0 samples the upper endpoint; 0 / 0 contributes nothing.
+  const double energy =
+      mass == 0 ? (kineticEnergy > 0 ? std::numeric_limits<double>::infinity()
+                                     : std::numeric_limits<double>::quiet_NaN())
+                : kineticEnergy / mass;
+  auto values = Lookup(charge, energy);
+  if (values.empty())
+    return;
+  const double edep =
+      step->GetTotalEnergyDeposit() / CLHEP::MeV * track->GetWeight();
+  ImageType4D::IndexType index4;
+  for (unsigned int i = 0; i < 3; ++i)
+    index4[i] = index[i];
+  // One actor-owned lock protects all shared accumulator updates.
+  std::lock_guard<std::mutex> lock(fScoringMutex);
+  ImageAddValue<ImageType>(cpp_amdm_restricted_edep_image, index, edep);
+  // LUT values are gamma[0..bins) followed by delta[0..bins).
+  // Keep both numerators raw; Python derives ratios only after worker scoring.
+  for (int b = 0; b < fBins; ++b) {
+    index4[3] = b;
+    const double delta = values[fBins + b] * edep;
+    ImageAddValue<ImageType4D>(cpp_amdm_delta_image, index4, delta);
+    ImageAddValue<ImageType4D>(cpp_amdm_gamma_image, index4, values[b] * delta);
+  }
+}

core/opengate_core/opengate_lib/GateAMDMActor.h

@@ -0,0 +1,61 @@
+/* --------------------------------------------------
+   Copyright (C): OpenGATE Collaboration
+   LGPL, see LICENSE.md.
+   Original AMDM implementation: Hermann Fuchs, April 2025.
+   -------------------------------------------------- */
+#ifndef GateAMDMActor_h
+#define GateAMDMActor_h
+
+#include "GateAMDMActor_lookUpTable.hh"
+#include "GateVActor.h"
+#include <G4RotationMatrix.hh>
+#include <itkImage.h>
+#include <mutex>
+
+/** Score weighted AMDM raw sums on a shared spatial grid with a separate bin
+ * axis. */
+class GateAMDMActor : public GateVActor {
+public:
+  using ImageType = itk::Image<double, 3>;
+  using ImageType4D = itk::Image<double, 4>;
+  /** Construct an MT-capable actor from the current Python user-info
+   * dictionary. */
+  explicit GateAMDMActor(py::dict &user_info);
+  /** Read voxel settings and validate/load the charge-energy LUT before
+   * scoring. */
+  void InitializeUserInfo(py::dict &user_info) override;
+  /** Initialize the base actor and create the 3D/4D ITK image holders. */
+  void InitializeCpp() override;
+  /** Reset the event count and attach fresh buffers to the current volume pose.
+   */
+  void BeginOfRunActionMasterThread(int run_id) override;
+  /** Count an event under the same actor-owned lock used for worker scoring. */
+  void BeginOfEventAction(const G4Event *) override;
+  /** Accumulate a step's weighted MeV, delta*MeV and gamma*delta*MeV raw sums.
+   */
+  void SteppingAction(G4Step *) override;
+  /** Return the physical-volume name used for run-time geometry attachment. */
+  std::string GetPhysicalVolumeName() const { return fPhysicalVolumeName; }
+  /** Select the resolved physical volume, including an inherited copy index. */
+  void SetPhysicalVolumeName(std::string name) { fPhysicalVolumeName = name; }
+  /** Return gamma bins then delta bins for charge and energy in MeV/n.
+   * Missing charges and NaN energies return an empty vector; endpoints clamp.
+   */
+  std::vector<double> Lookup(int charge, double energy) const;
+
+  ImageType::Pointer cpp_amdm_restricted_edep_image;
+  ImageType4D::Pointer
+      cpp_amdm_delta_image; // raw sums, never normalized in place
+  ImageType4D::Pointer cpp_amdm_gamma_image;
+  int NbOfEvent = 0;
+
+private:
+  GateAMDMLookUpTable fTable;
+  int fBins = 10;
+  std::string fHitType;
+  std::string fPhysicalVolumeName;
+  G4ThreeVector fTranslation;
+  G4RotationMatrix fRotation;
+  std::mutex fScoringMutex; // owned by this actor, shared by its workers
+};
+#endif

core/opengate_core/opengate_lib/GateAMDMActor_lookUpTable.cpp

@@ -0,0 +1,100 @@
+/* Copyright (C): OpenGATE Collaboration. See LICENSE.md for LGPL terms. */
+#include "GateAMDMActor_lookUpTable.hh"
+#include <algorithm>
+#include <cmath>
+#include <fstream>
+#include <limits>
+#include <sstream>
+#include <stdexcept>
+
+void GateAMDMLookUpTable::Read(const std::string &filename, int bins) {
+  if (bins <= 0 || bins > (std::numeric_limits<int>::max() - 2) / 2)
+    throw std::invalid_argument("AMDM_Bins must be a positive integer");
+  std::ifstream input(filename);
+  if (!input)
+    throw std::invalid_argument("AMDM cannot open LUT: " + filename);
+  // Parse transactionally: failed reinitialization must not leave a partial
+  // LUT.
+  std::map<int, std::vector<Row>> data;
+  std::string line;
+  int previousCharge = 0;
+  size_t lineNumber = 0;
+  while (std::getline(input, line)) {
+    ++lineNumber;
+    line = line.substr(0, line.find('#'));
+    std::istringstream stream(line);
+    std::vector<double> columns;
+    std::string token;
+    // Include the original input line in every validation error.
+    auto fail = [&](const std::string &reason) {
+      throw std::invalid_argument("AMDM LUT " + filename + ":" +
+                                  std::to_string(lineNumber) + ": " + reason);
+    };
+    while (stream >> token) {
+      try {
+        size_t used;
+        double value = std::stod(token, &used);
+        if (used != token.size() || !std::isfinite(value))
+          fail("entries must be finite numeric values");
+        columns.push_back(value);
+      } catch (const std::exception &) {
+        fail("entries must be finite numeric values");
+      }
+    }
+    if (columns.empty())
+      continue;
+    if (columns.size() != static_cast<size_t>(2 + 2 * bins))
+      fail("column count must equal 2 + 2 * AMDM_Bins");
+    double key = columns[0];
+    if (key <= 0 || key > std::numeric_limits<int>::max() ||
+        std::floor(key) != key)
+      fail("charge must be a positive integer");
+    int charge = static_cast<int>(key);
+    if (charge < previousCharge)
+      fail("charge groups must be increasing");
+    if (columns[1] <= 0)
+      fail("energy keys must be positive (MeV/n)");
+    auto &group = data[charge];
+    if (!group.empty() && columns[1] <= group.back().energy)
+      fail(
+          "energy keys must be unique and strictly increasing within a charge");
+    group.push_back({columns[1], {columns.begin() + 2, columns.end()}});
+    previousCharge = charge;
+  }
+  if (input.bad())
+    throw std::invalid_argument("AMDM error reading LUT: " + filename);
+  if (data.empty())
+    throw std::invalid_argument("AMDM LUT contains no data: " + filename);
+  fData = std::move(data);
+}
+
+bool GateAMDMLookUpTable::Find(int charge, double energy,
+                               std::vector<double> &values) const {
+  values.clear();
+  auto found = fData.find(charge);
+  if (found == fData.end() || std::isnan(energy))
+    return false;
+  const auto &rows = found->second;
+  // Clamp before inspecting a lower_bound iterator, including the final charge.
+  if (energy <= rows.front().energy) {
+    values = rows.front().values;
+  } else if (energy >= rows.back().energy) {
+    values = rows.back().values;
+  } else {
+    auto upper = std::lower_bound(
+        rows.begin(), rows.end(), energy,
+        [](const Row &row, double value) { return row.energy < value; });
+    if (upper->energy == energy) {
+      values = upper->values;
+    } else {
+      const auto &lower = *std::prev(upper);
+      const double fraction =
+          (energy - lower.energy) / (upper->energy - lower.energy);
+      values.resize(lower.values.size());
+      for (size_t b = 0; b < values.size(); ++b)
+        values[b] =
+            lower.values[b] + fraction * (upper->values[b] - lower.values[b]);
+    }
+  }
+  return true;
+}

core/opengate_core/opengate_lib/GateAMDMActor_lookUpTable.hh

@@ -0,0 +1,31 @@
+/* Copyright (C): OpenGATE Collaboration. See LICENSE.md for LGPL terms. */
+#ifndef GateAMDMActor_lookUpTable_h
+#define GateAMDMActor_lookUpTable_h
+
+#include <map>
+#include <string>
+#include <vector>
+
+/** Immutable during scoring; each charge owns its energy grid and endpoints. */
+class GateAMDMLookUpTable {
+public:
+  /** Load a whitespace/comment table with 2 + 2*bins finite columns per row.
+   * Validate positive keys and ordered unique energy grids. Raise
+   * std::invalid_argument on malformed/unreadable input and retain old data.
+   */
+  void Read(const std::string &filename, int bins);
+  /** Fill values with gamma bins followed by delta bins at energy in MeV/n.
+   * Interpolate linearly within a group and clamp at its endpoints. Return
+   * false and clear values for a missing charge or NaN energy; infinities
+   * clamp.
+   */
+  bool Find(int charge, double energy, std::vector<double> &values) const;
+
+private:
+  struct Row {
+    double energy;
+    std::vector<double> values; // gamma bins, then delta bins
+  };
+  std::map<int, std::vector<Row>> fData;
+};
+#endif