feat: implemented a wrapper around the UserSpecialCuts for EkinMin but with additional filtering for particle type

Author: MoritzNeuberger

include/RMGHardware.hh

@@ -20,6 +20,7 @@
 #include <memory>
 #include <set>
 #include <string>
+#include <unordered_map>
 #include <vector>
 
 #include "G4GenericMessenger.hh"
@@ -33,6 +34,7 @@
 #include "RMGVOutputScheme.hh"
 
 /** @brief Class to handle the detector geometry hardware, extends @c G4VUserDetectorConstruction. */
+class G4LogicalVolume;
 class G4VPhysicalVolume;
 class RMGHardware : public G4VUserDetectorConstruction {
 
@@ -158,6 +160,23 @@ class RMGHardware : public G4VUserDetectorConstruction {
      */
     void SetMaxStepLimit(double max_step, std::string name);
 
+    /** @brief Set a minimum kinetic energy cut for one selected particle in a volume.
+     *
+     * @details This sets @c G4UserLimits::SetUserMinEkine on matching logical volumes, but the
+     * corresponding cut process is applied only to registered particles.
+     *
+     * @param ekin_min Minimum kinetic energy.
+     * @param name Physical volume name or regex.
+     * @param particle_name Geant4 particle name.
+     */
+    void SetEminLimitForParticle(double ekin_min, std::string name, std::string particle_name);
+
+    /** @brief Check if selective EkinMin should be applied for the current track context. */
+    static bool IsEminLimitParticleSelected(
+        const G4LogicalVolume* logical,
+        const std::string& particle_name
+    );
+
     void PrintListOfLogicalVolumes() { RMGNavigationTools::PrintListOfLogicalVolumes(); }
     void PrintListOfPhysicalVolumes() { RMGNavigationTools::PrintListOfPhysicalVolumes(); }
 
@@ -176,6 +195,13 @@ class RMGHardware : public G4VUserDetectorConstruction {
     /// Mapping between physical volume names and maximum (user) step size to apply
     std::map<std::string, double> fPhysVolStepLimits;
 
+    struct RMGSelectiveEminLimit {
+        double ekin_min;
+        std::set<std::string> particles;
+    };
+    std::map<std::string, RMGSelectiveEminLimit> fPhysVolEminLimits;
+    static std::unordered_map<const G4LogicalVolume*, std::set<std::string>> fLogicalVolEminParticles;
+
     // one element for each sensitive detector physical volume
     std::map<std::pair<std::string, int>, RMGDetectorMetadata> fDetectorMetadata;
 

include/RMGHardwareMessenger.hh

@@ -45,14 +45,19 @@ class RMGHardwareMessenger : public G4UImessenger {
     /** @brief Define the commands to set the step limits. */
     void DefineStepLimits();
 
+    /** @brief Define the commands to set particle-selective minimum kinetic energy limits. */
+    void DefineSelectiveEminLimits();
+
   private:
 
     RMGHardware* fHardware;
     G4UIcommand* fRegisterCmd = nullptr;
     G4UIcmdWithADoubleAndUnit* fStepLimitsCmd = nullptr;
+    G4UIcommand* fSelectiveEminLimitCmd = nullptr;
 
     void RegisterDetectorCmd(const std::string& parameters);
     void StepLimitsCmd(const std::string& parameters);
+    void SelectiveEminLimitCmd(const std::string& parameters);
 };
 
 #endif

include/RMGSelectiveEkinMinCutProcess.hh

@@ -0,0 +1,27 @@
+// Copyright (C) 2026
+//
+// This program is free software: you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any
+// later version.
+
+#ifndef _RMG_SELECTIVE_EKIN_MIN_CUT_PROCESS_HH_
+#define _RMG_SELECTIVE_EKIN_MIN_CUT_PROCESS_HH_
+
+#include "G4UserSpecialCuts.hh"
+
+class RMGSelectiveEkinMinCutProcess : public G4UserSpecialCuts {
+
+  public:
+
+    explicit RMGSelectiveEkinMinCutProcess(const G4String& process_name = "UserSpecialCut");
+    ~RMGSelectiveEkinMinCutProcess() override = default;
+
+    G4double PostStepGetPhysicalInteractionLength(
+        const G4Track& aTrack,
+        G4double previousStepSize,
+        G4ForceCondition* condition
+    ) override;
+};
+
+#endif

src/CMakeLists.txt

@@ -41,6 +41,7 @@ set(PROJECT_PUBLIC_HEADERS
     ${_root}/include/RMGRun.hh
     ${_root}/include/RMGScintillatorDetector.hh
     ${_root}/include/RMGScintillatorOutputScheme.hh
+    ${_root}/include/RMGSelectiveEkinMinCutProcess.hh
     ${_root}/include/RMGStackingAction.hh
     ${_root}/include/RMGSteppingAction.hh
     ${_root}/include/RMGTools.hh
@@ -93,6 +94,7 @@ set(PROJECT_SOURCES
     ${_root}/src/RMGRunAction.cc
     ${_root}/src/RMGScintillatorDetector.cc
     ${_root}/src/RMGScintillatorOutputScheme.cc
+    ${_root}/src/RMGSelectiveEkinMinCutProcess.cc
     ${_root}/src/RMGStackingAction.cc
     ${_root}/src/RMGSteppingAction.cc
     ${_root}/src/RMGTrackingAction.cc

src/RMGHardware.cc

@@ -15,6 +15,7 @@
 
 #include "RMGHardware.hh"
 
+#include <algorithm>
 #include <filesystem>
 #include <format>
 namespace fs = std::filesystem;
@@ -42,6 +43,8 @@ namespace fs = std::filesystem;
 #include "RMGTools.hh"
 #include "RMGVertexOutputScheme.hh"
 
+#include "fmt/ranges.h"
+
 #if RMG_HAS_GDML
 #include "G4GDMLParser.hh"
 #endif
@@ -52,6 +55,9 @@ G4ThreadLocal std::vector<std::shared_ptr<RMGVOutputScheme>> RMGHardware::fActiv
 
 G4ThreadLocal bool RMGHardware::fActiveDetectorsInitialized = false;
 
+std::unordered_map<const G4LogicalVolume*, std::set<std::string>>
+    RMGHardware::fLogicalVolEminParticles = {};
+
 RMGHardware::RMGHardware() { this->DefineCommands(); }
 
 G4VPhysicalVolume* RMGHardware::Construct() {
@@ -166,11 +172,65 @@ G4VPhysicalVolume* RMGHardware::Construct() {
       );
     } else {
       for (const auto& vol : volumes) {
-        vol->GetLogicalVolume()->SetUserLimits(new G4UserLimits(el.second));
+        G4LogicalVolume* logical = vol->GetLogicalVolume();
+        auto userLimits = logical->GetUserLimits();
+        if (!userLimits) { userLimits = new G4UserLimits(); }
+        userLimits->SetMaxAllowedStep(el.second);
+        logical->SetUserLimits(userLimits);
       }
     }
   }
 
+  // attach particle-selective user min kinetic energy cuts to logical volumes
+  fLogicalVolEminParticles.clear();
+  std::unordered_map<const G4LogicalVolume*, double> logical_ekin_limits;
+  for (const auto& [vol_name, cfg] : fPhysVolEminLimits) {
+    RMGLog::OutFormat(
+        RMGLog::debug,
+        "Setting selective min user kinetic energy for volume '{}' to {} ({})",
+        vol_name,
+        cfg.ekin_min,
+        fmt::join(cfg.particles, ",")
+    );
+    auto volumes = RMGNavigationTools::FindPhysicalVolume(vol_name);
+    if (volumes.empty()) {
+      RMGLog::Out(
+          RMGLog::error,
+          "No matching volumes for '{}' found, skipping selective user ekin min limit setting",
+          vol_name
+      );
+      continue;
+    }
+
+    for (const auto& vol : volumes) {
+      G4LogicalVolume* logical = vol->GetLogicalVolume();
+
+      const auto ekin_it = logical_ekin_limits.find(logical);
+      if (ekin_it != logical_ekin_limits.end() && ekin_it->second != cfg.ekin_min) {
+        RMGLog::OutFormat(
+            RMGLog::error,
+            "Conflicting selective ekin min limits for logical volume {} ({} vs {} MeV), "
+            "skipping registration from pattern '{}'",
+            logical->GetName(),
+            ekin_it->second,
+            cfg.ekin_min,
+            vol_name
+        );
+        continue;
+      }
+
+      auto userLimits = logical->GetUserLimits();
+      if (!userLimits) userLimits = new G4UserLimits();
+      userLimits->SetUserMinEkine(cfg.ekin_min);
+      logical->SetUserLimits(userLimits);
+
+      logical_ekin_limits.emplace(logical, cfg.ekin_min);
+
+      auto [it, inserted] = fLogicalVolEminParticles.try_emplace(logical, cfg.particles);
+      if (!inserted) { it->second.insert(cfg.particles.begin(), cfg.particles.end()); }
+    }
+  }
+
   // register staged detectors now
   if (!fStagedDetectors.empty()) {
     RMGLog::Out(RMGLog::debug, "Registering staged detectors");
@@ -489,6 +549,48 @@ void RMGHardware::SetMaxStepLimit(double max_step, std::string name) {
   RMGLog::OutFormat(RMGLog::detail, "Set step limits for {:s} to {:.2f} mm", name, max_step);
 }
 
+void RMGHardware::SetEminLimitForParticle(double ekin_min, std::string name, std::string particle_name) {
+  if (particle_name.empty()) {
+    RMGLog::Out(
+        RMGLog::error,
+        "No particle name was provided for selective ekin min limit, ignoring command"
+    );
+    return;
+  }
+
+  auto [it, inserted] = fPhysVolEminLimits.try_emplace(name, RMGSelectiveEminLimit{ekin_min, {}});
+  if (!inserted && it->second.ekin_min != ekin_min) {
+    RMGLog::OutFormat(
+        RMGLog::error,
+        "Conflicting selective ekin min limit for volume {:s}: existing {:.2f} MeV, requested "
+        "{:.2f} MeV",
+        name,
+        it->second.ekin_min,
+        ekin_min
+    );
+    return;
+  }
+
+  it->second.particles.insert(particle_name);
+  RMGLog::OutFormat(
+      RMGLog::detail,
+      "Set selective ekin min limit for {:s} to {:.2f} MeV (particle: {})",
+      name,
+      ekin_min,
+      particle_name
+  );
+}
+
+bool RMGHardware::IsEminLimitParticleSelected(
+    const G4LogicalVolume* logical,
+    const std::string& particle_name
+) {
+  if (!logical) return false;
+  const auto it = fLogicalVolEminParticles.find(logical);
+  if (it == fLogicalVolEminParticles.end()) return false;
+  return it->second.find(particle_name) != it->second.end();
+}
+
 void RMGHardware::RegisterDetectorsFromGDML(std::string s) {
   if (s == "All") {
     for (const auto dt : magic_enum::enum_values<RMGDetectorType>()) {

src/RMGHardwareMessenger.cc

@@ -25,14 +25,20 @@ RMGHardwareMessenger::RMGHardwareMessenger(RMGHardware* hw) : fHardware(hw) {
 
   this->DefineRegisterDetector();
   this->DefineStepLimits();
+  this->DefineSelectiveEminLimits();
 }
 
-RMGHardwareMessenger::~RMGHardwareMessenger() { delete fRegisterCmd; }
+RMGHardwareMessenger::~RMGHardwareMessenger() {
+  delete fRegisterCmd;
+  delete fStepLimitsCmd;
+  delete fSelectiveEminLimitCmd;
+}
 
 void RMGHardwareMessenger::SetNewValue(G4UIcommand* command, G4String newValues) {
 
   if (command == fRegisterCmd) RegisterDetectorCmd(newValues);
   else if (command == fStepLimitsCmd) StepLimitsCmd(newValues);
+  else if (command == fSelectiveEminLimitCmd) SelectiveEminLimitCmd(newValues);
 }
 
 void RMGHardwareMessenger::DefineRegisterDetector() {
@@ -89,6 +95,33 @@ void RMGHardwareMessenger::DefineStepLimits() {
   fStepLimitsCmd->AvailableForStates(G4State_PreInit);
 }
 
+void RMGHardwareMessenger::DefineSelectiveEminLimits() {
+
+  fSelectiveEminLimitCmd = new G4UIcommand("/RMG/Geometry/SetEkinMinForParticle", this);
+  fSelectiveEminLimitCmd->SetGuidance(
+      "Sets minimum kinetic energy for one selected particle in a detector volume"
+  );
+
+  auto p_ekin = new G4UIparameter("ekin_min", 'd', false);
+  p_ekin->SetGuidance("minimum kinetic energy");
+  fSelectiveEminLimitCmd->SetParameter(p_ekin);
+
+  auto p_unit = new G4UIparameter("unit", 's', false);
+  p_unit->SetGuidance("energy unit");
+  p_unit->SetParameterCandidates("eV keV MeV GeV TeV");
+  fSelectiveEminLimitCmd->SetParameter(p_unit);
+
+  auto p_pv = new G4UIparameter("pv_name", 's', false);
+  p_pv->SetGuidance("Detector physical volume, accepts regex patterns");
+  fSelectiveEminLimitCmd->SetParameter(p_pv);
+
+  auto p_particle = new G4UIparameter("particle_name", 's', false);
+  p_particle->SetGuidance("Geant4 particle name, e.g. e-, e+, gamma");
+  fSelectiveEminLimitCmd->SetParameter(p_particle);
+
+  fSelectiveEminLimitCmd->AvailableForStates(G4State_PreInit);
+}
+
 void RMGHardwareMessenger::RegisterDetectorCmd(const std::string& parameters) {
   G4Tokenizer next(parameters);
 
@@ -117,4 +150,16 @@ void RMGHardwareMessenger::StepLimitsCmd(const std::string& parameters) {
   fHardware->SetMaxStepLimit(num, pv_name);
 }
 
+void RMGHardwareMessenger::SelectiveEminLimitCmd(const std::string& parameters) {
+  G4Tokenizer next(parameters);
+
+  auto x = next();
+  auto y = next();
+  auto num = G4UIcmdWithADoubleAndUnit::GetNewDoubleValue((x + " " + y).c_str());
+  auto pv_name = next();
+  auto particle_name = next();
+
+  fHardware->SetEminLimitForParticle(num, pv_name, particle_name);
+}
+
 // vim: tabstop=2 shiftwidth=2 expandtab

src/RMGPhysics.cc

@@ -70,13 +70,15 @@
 #include "G4StepLimiterPhysics.hh"
 #include "G4StoppingPhysics.hh"
 #include "G4ThermalNeutrons.hh"
+#include "G4UserSpecialCuts.hh"
 #include "G4Version.hh"
 
 #include "RMGConfig.hh"
 #include "RMGInnerBremsstrahlungProcess.hh"
 #include "RMGLog.hh"
 #include "RMGNeutronCaptureProcess.hh"
 #include "RMGOpWLSProcess.hh"
+#include "RMGSelectiveEkinMinCutProcess.hh"
 #include "RMGTools.hh"
 
 namespace u = CLHEP;
@@ -360,6 +362,21 @@ void RMGPhysics::ConstructProcess() {
   // add step limits
   auto step_limits = new G4StepLimiterPhysics();
   step_limits->ConstructProcess();
+
+  // replace generic UserSpecialCuts with particle-selective variant used by
+  // /RMG/Geometry/SetEkinMinForParticle.
+  GetParticleIterator()->reset();
+  while ((*GetParticleIterator())()) {
+    auto particle = GetParticleIterator()->value();
+    auto proc_manager = particle->GetProcessManager();
+    if (!proc_manager) continue;
+
+    auto user_special_cut = proc_manager->GetProcess("UserSpecialCut");
+    if (!user_special_cut) continue;
+
+    proc_manager->RemoveProcess(user_special_cut);
+    proc_manager->AddDiscreteProcess(new RMGSelectiveEkinMinCutProcess());
+  }
 }
 
 void RMGPhysics::ConstructOptical() {