Merge pull request #121 from ManuelHu/isotope-filter

add optional isotope filter output scheme

docs/g4manual2rst.py

@@ -35,6 +35,7 @@ def remove_whitespace_lines_end(lines: list):
 
 idx = 0
 in_cmdblock = False
+in_guidance = False
 lastlevel = -1
 leveldiff = 0
 
@@ -52,16 +53,23 @@ def remove_whitespace_lines_end(lines: list):
         lastlevel = -1
         leveldiff = 0
     elif in_cmdblock and (line == "Guidance :"):
-        pass
-    elif in_cmdblock and (inlines[idx - 1] == "Guidance :") and not line.startswith(" " * 2):
-        outlines.extend([line, ""])
+        in_guidance = True
+    elif in_cmdblock and in_guidance and not line.startswith(" "):
+        if line.startswith("note: "):
+            outlines.extend([".. note ::", "", (" " * 4) + line[6:], ""])
+        elif line.strip() == "":
+            in_guidance = False
+        elif line != "":
+            outlines.extend([line, ""])
     elif in_cmdblock and line == " Commands : " and not inlines[idx + 1].startswith(" " * 3):
         # ignore directories with no commands.
         pass
     elif in_cmdblock and line == " Sub-directories : " and inlines[idx + 1] == " Commands : ":
         # ignore directories with no commands.
         pass
     elif in_cmdblock and line != "":
+        in_guidance = False
+
         stripped_line = line.lstrip()
         indent = math.ceil((len(line) - len(stripped_line)) / 2)
         if line.startswith(" Range of parameters :"):

include/RMGIsotopeFilterOutputScheme.hh

@@ -0,0 +1,62 @@
+// Copyright (C) 2022 Luigi Pertoldi <[email protected]>
+//
+// 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.
+//
+// This program is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+#ifndef _RMG_ISOTOPE_FILTER_OUTPUT_SCHEME_HH_
+#define _RMG_ISOTOPE_FILTER_OUTPUT_SCHEME_HH_
+
+#include <optional>
+#include <set>
+
+#include "G4AnalysisManager.hh"
+#include "G4GenericMessenger.hh"
+#include "G4VUserEventInformation.hh"
+
+#include "RMGVOutputScheme.hh"
+
+class RMGIsotopeFilterEventInformation : public G4VUserEventInformation {
+
+  public:
+
+    RMGIsotopeFilterEventInformation() = default;
+    inline void Print() const override {}
+};
+
+class G4Event;
+class RMGIsotopeFilterOutputScheme : public RMGVOutputScheme {
+
+  public:
+
+    RMGIsotopeFilterOutputScheme();
+
+    bool ShouldDiscardEvent(const G4Event*) override;
+    std::optional<bool> StackingActionNewStage(int) override;
+    std::optional<G4ClassificationOfNewTrack> StackingActionClassify(const G4Track*, int) override;
+    void TrackingActionPre(const G4Track* aTrack) override;
+
+    inline void AddIsotope(int a, int z) { fIsotopes.insert({a, z}); }
+
+  private:
+
+    std::unique_ptr<G4GenericMessenger> fMessenger;
+    void DefineCommands();
+
+    std::set<std::pair<int, int>> fIsotopes;
+
+    bool fDiscardPhotonsIfIsotopeNotProduced = false;
+};
+
+#endif
+
+// vim: tabstop=2 shiftwidth=2 expandtab

include/RMGUserInit.hh

@@ -26,6 +26,7 @@
 #include "G4UserSteppingAction.hh"
 #include "G4UserTrackingAction.hh"
 
+#include "RMGIsotopeFilterOutputScheme.hh"
 #include "RMGLog.hh"
 #include "RMGVOutputScheme.hh"
 
@@ -65,6 +66,11 @@ class RMGUserInit {
       Add<T>(&fOptionalOutputSchemes, name, std::forward<Args>(args)...);
     }
 
+    // default output schemes
+    inline void RegisterDefaultOptionalOutputSchemes() {
+      AddOptionalOutputScheme<RMGIsotopeFilterOutputScheme>("IsotopeFilter");
+    }
+
     inline void ActivateOptionalOutputScheme(std::string name) {
       auto it = fOptionalOutputSchemes.find(name);
       if (it == fOptionalOutputSchemes.end()) {

include/RMGVOutputScheme.hh

@@ -49,7 +49,7 @@ class RMGVOutputScheme {
     virtual inline std::optional<bool> StackingActionNewStage(const int) { return std::nullopt; }
 
     // hook into G4TrackingAction
-    virtual inline void TrackingActionPre(const G4Track* aTrack) {};
+    virtual inline void TrackingActionPre(const G4Track*) {};
 
     inline void SetNtuplePerDetector(bool ntuple_per_det) { fNtuplePerDetector = ntuple_per_det; }
 

src/CMakeLists.txt

@@ -13,6 +13,7 @@ set(PROJECT_PUBLIC_HEADERS
     ${_root}/include/RMGGeneratorUtil.hh
     ${_root}/include/RMGGermaniumDetector.hh
     ${_root}/include/RMGGermaniumOutputScheme.hh
+    ${_root}/include/RMGIsotopeFilterOutputScheme.hh
     ${_root}/include/RMGLog.hh
     ${_root}/include/RMGManager.hh
     ${_root}/include/RMGMasterGenerator.hh
@@ -48,6 +49,7 @@ set(PROJECT_SOURCES
     ${_root}/src/RMGGeneratorUtil.cc
     ${_root}/src/RMGGermaniumDetector.cc
     ${_root}/src/RMGGermaniumOutputScheme.cc
+    ${_root}/src/RMGIsotopeFilterOutputScheme.cc
     ${_root}/src/RMGLog.cc
     ${_root}/src/RMGManager.cc
     ${_root}/src/RMGMasterGenerator.cc

src/RMGGermaniumOutputScheme.cc

@@ -19,6 +19,7 @@
 
 #include "G4AnalysisManager.hh"
 #include "G4Event.hh"
+#include "G4EventManager.hh"
 #include "G4HCtable.hh"
 #include "G4OpticalPhoton.hh"
 #include "G4SDManager.hh"
@@ -163,19 +164,24 @@ void RMGGermaniumOutputScheme::StoreEvent(const G4Event* event) {
 
 std::optional<G4ClassificationOfNewTrack> RMGGermaniumOutputScheme::StackingActionClassify(const G4Track* aTrack,
     int stage) {
+  // we are only interested in stacking optical photons into stage 1 after stage 0 finished.
   if (stage != 0) return std::nullopt;
-  // defer tracking of optical photons.
+
+  // defer tracking of optical photons, irrespective of our settings.
   if (aTrack->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) return fWaiting;
-  return fUrgent;
+  return std::nullopt;
 }
 
 std::optional<bool> RMGGermaniumOutputScheme::StackingActionNewStage(const int stage) {
+  // we are only interested in stacking optical photons into stage 1 after stage 0 finished.
   if (stage != 0) return std::nullopt;
-  if (!fDiscardPhotonsIfNoGermaniumEdep) return true;
-  auto run_man = RMGManager::Instance()->GetG4RunManager();
-  const auto event = run_man->GetCurrentEvent();
-  // discard all waiting events, as there was no energy deposition in Germanium.
-  return !ShouldDiscardEvent(event);
+  // if we do not want to discard any photons ourselves, let other output schemes decide (i.e. not
+  // force `true` on them).
+  if (!fDiscardPhotonsIfNoGermaniumEdep) return std::nullopt;
+
+  const auto event = G4EventManager::GetEventManager()->GetConstCurrentEvent();
+  // discard all waiting events, if there was no energy deposition in Germanium.
+  return ShouldDiscardEvent(event) ? std::make_optional(false) : std::nullopt;
 }
 
 void RMGGermaniumOutputScheme::DefineCommands() {
@@ -201,8 +207,10 @@ void RMGGermaniumOutputScheme::DefineCommands() {
       .SetStates(G4State_Idle);
 
   fMessenger->DeclareProperty("DiscardPhotonsIfNoGermaniumEdep", fDiscardPhotonsIfNoGermaniumEdep)
-      .SetGuidance(
-          "Discard optical photons (before simulating them), if no edep in germanium detectors.")
+      .SetGuidance("Discard optical photons (before simulating them), if no edep in germanium "
+                   "detectors occurred in the same event.")
+      .SetGuidance("note: If another output scheme also requests the photons to be discarded, the "
+                   "germanium edep filter does not force the photons to be simulated.")
       .SetStates(G4State_Idle);
 }
 

src/RMGIsotopeFilterOutputScheme.cc

@@ -0,0 +1,105 @@
+// Copyright (C) 2022 Luigi Pertoldi <[email protected]>
+//
+// 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.
+//
+// This program is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+#include "RMGIsotopeFilterOutputScheme.hh"
+
+#include <set>
+
+#include "G4Event.hh"
+#include "G4EventManager.hh"
+#include "G4OpticalPhoton.hh"
+
+#include "RMGLog.hh"
+
+RMGIsotopeFilterOutputScheme::RMGIsotopeFilterOutputScheme() { this->DefineCommands(); }
+
+void RMGIsotopeFilterOutputScheme::TrackingActionPre(const G4Track* aTrack) {
+  const auto particle = aTrack->GetParticleDefinition();
+  if (!particle->IsGeneralIon()) return;
+  const int z = particle->GetAtomicNumber();
+  const int a = particle->GetAtomicMass();
+  if (z == 0 || a == 0) return; // not an isotope, but any other particle.
+
+  if (fIsotopes.find({a, z}) == fIsotopes.end()) return;
+
+  auto event = G4EventManager::GetEventManager()->GetNonconstCurrentEvent();
+  auto info = event->GetUserInformation();
+  if (info != nullptr && dynamic_cast<RMGIsotopeFilterEventInformation*>(info) == nullptr) {
+    RMGLog::Out(RMGLog::error,
+        "other user info class found, instead of RMGIsotopeFilterEventInformation");
+    return;
+  }
+  if (info == nullptr) event->SetUserInformation(new RMGIsotopeFilterEventInformation());
+}
+
+// invoked in RMGEventAction::EndOfEventAction()
+bool RMGIsotopeFilterOutputScheme::ShouldDiscardEvent(const G4Event* event) {
+  // exit fast if no threshold is configured.
+  if (fIsotopes.empty()) return false;
+
+  auto info = event->GetUserInformation();
+  if (info != nullptr && dynamic_cast<RMGIsotopeFilterEventInformation*>(info) == nullptr) {
+    // Someone else tries to set Event user information.
+    RMGLog::Out(RMGLog::error,
+        "other user info class found, instead of RMGIsotopeFilterEventInformation");
+    return false;
+  }
+  return info == nullptr;
+}
+
+std::optional<G4ClassificationOfNewTrack> RMGIsotopeFilterOutputScheme::
+    StackingActionClassify(const G4Track* aTrack, int stage) {
+  // we are only interested in stacking optical photons into stage 1 after stage 0 finished.
+  if (stage != 0) return std::nullopt;
+
+  // defer tracking of optical photons.
+  if (aTrack->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) return fWaiting;
+  return std::nullopt;
+}
+
+std::optional<bool> RMGIsotopeFilterOutputScheme::StackingActionNewStage(const int stage) {
+  // we are only interested in stacking optical photons into stage 1 after stage 0 finished.
+  if (stage != 0) return std::nullopt;
+  // if we do not want to discard any photons ourselves, let other output schemes decide (i.e. not
+  // force `true` on them).
+  if (!fDiscardPhotonsIfIsotopeNotProduced) return std::nullopt;
+
+  const auto event = G4EventManager::GetEventManager()->GetConstCurrentEvent();
+  // discard all waiting events, if there were none of the requested isotopes produced.
+  return ShouldDiscardEvent(event) ? std::make_optional(false) : std::nullopt;
+}
+
+void RMGIsotopeFilterOutputScheme::DefineCommands() {
+
+  fMessenger = std::make_unique<G4GenericMessenger>(this, "/RMG/Output/IsotopeFilter/",
+      "Commands for filtering event out by created isotopes.");
+
+  fMessenger->DeclareMethod("AddIsotope", &RMGIsotopeFilterOutputScheme::AddIsotope)
+      .SetGuidance("Add an isotope to the list. Only events that have a track with this isotope at "
+                   "any point in time will be persisted.")
+      .SetParameterName(0, "A", false, false)
+      .SetParameterName(1, "Z", false, false)
+      .SetStates(G4State_Idle);
+
+  fMessenger
+      ->DeclareProperty("DiscardPhotonsIfIsotopeNotProduced", fDiscardPhotonsIfIsotopeNotProduced)
+      .SetGuidance("Discard optical photons (before simulating them), if the specified isotopes "
+                   "had not been produced in the same event.")
+      .SetGuidance("note: If another output scheme also requests the photons to be discarded, the "
+                   "isotope filter does not force the photons to be simulated.")
+      .SetStates(G4State_Idle);
+}
+
+// vim: tabstop=2 shiftwidth=2 expandtab

src/RMGManager.cc

@@ -66,7 +66,9 @@ RMGManager::RMGManager(std::string app_name, int argc, char** argv)
   // limit Geant4 stacktrace dumping to segfaults
   G4Backtrace::DefaultSignals() = std::set<int>{SIGSEGV};
 
+  // initialize central hook for dependent applications to influence our output.
   fUserInit = std::make_shared<RMGUserInit>();
+  fUserInit->RegisterDefaultOptionalOutputSchemes();
 
   this->DefineCommands();
 }