add option to discard event with energy threshold in germanium

include/RMGGermaniumDetector.hh

@@ -77,9 +77,6 @@ class RMGGermaniumDetector : public G4VSensitiveDetector {
   private:
 
     RMGGermaniumDetectorHitsCollection* fHitsCollection = nullptr;
-
-    std::unique_ptr<G4GenericMessenger> fMessenger = nullptr;
-    void DefineCommands();
 };
 
 extern G4ThreadLocal G4Allocator<RMGGermaniumDetectorHit>* RMGGermaniumDetectorHitAllocator;

include/RMGGermaniumOutputScheme.hh

@@ -16,17 +16,39 @@
 #ifndef _RMG_GERMANIUM_OUTPUT_SCHEME_HH_
 #define _RMG_GERMANIUM_OUTPUT_SCHEME_HH_
 
+#include <set>
+
+#include "G4AnalysisManager.hh"
+#include "G4GenericMessenger.hh"
+
+#include "RMGGermaniumDetector.hh"
 #include "RMGVOutputScheme.hh"
 
 class G4Event;
 class RMGGermaniumOutputScheme : public RMGVOutputScheme {
 
   public:
 
-    inline RMGGermaniumOutputScheme(G4AnalysisManager* ana_man) : RMGVOutputScheme(ana_man) {}
+    RMGGermaniumOutputScheme();
 
     void AssignOutputNames(G4AnalysisManager* ana_man) override;
-    void EndOfEventAction(const G4Event*) override;
+    void StoreEvent(const G4Event*) override;
+    bool ShouldDiscardEvent(const G4Event*) override;
+
+    inline void SetEdepCutLow(double threshold) { fEdepCutLow = threshold; }
+    inline void SetEdepCutHigh(double threshold) { fEdepCutHigh = threshold; }
+    inline void AddEdepCutDetector(int det_uid) { fEdepCutDetectors.insert(det_uid); }
+
+  private:
+
+    RMGGermaniumDetectorHitsCollection* GetHitColl(const G4Event*);
+
+    std::unique_ptr<G4GenericMessenger> fMessenger;
+    void DefineCommands();
+
+    double fEdepCutLow = -1;
+    double fEdepCutHigh = -1;
+    std::set<int> fEdepCutDetectors;
 };
 
 #endif

include/RMGHardware.hh

@@ -28,6 +28,7 @@
 
 #include "RMGHardwareMessenger.hh"
 #include "RMGNavigationTools.hh"
+#include "RMGVOutputScheme.hh"
 
 class G4VPhysicalVolume;
 class RMGHardware : public G4VUserDetectorConstruction {
@@ -63,6 +64,7 @@ class RMGHardware : public G4VUserDetectorConstruction {
       return fDetectorMetadata.at(det);
     }
     inline const auto& GetActiveDetectorList() { return fActiveDetectors; }
+    inline auto GetActiveOutputScheme(DetectorType type) { return fActiveOutputSchemes.at(type); }
 
     inline void IncludeGDMLFile(std::string filename) { fGDMLFiles.emplace_back(filename); }
     inline virtual G4VPhysicalVolume* DefineGeometry() { return nullptr; }
@@ -84,6 +86,7 @@ class RMGHardware : public G4VUserDetectorConstruction {
     // one element for each sensitive detector physical volume
     std::map<std::pair<std::string, int>, DetectorMetadata> fDetectorMetadata;
     std::set<DetectorType> fActiveDetectors;
+    std::map<DetectorType, std::shared_ptr<RMGVOutputScheme>> fActiveOutputSchemes;
     bool fActiveDetectorsInitialized = false;
 
     std::unique_ptr<G4GenericMessenger> fMessenger;

include/RMGOpticalDetector.hh

@@ -76,9 +76,6 @@ class RMGOpticalDetector : public G4VSensitiveDetector {
   private:
 
     RMGOpticalDetectorHitsCollection* fHitsCollection = nullptr;
-
-    std::unique_ptr<G4GenericMessenger> fMessenger = nullptr;
-    void DefineCommands();
 };
 
 extern G4ThreadLocal G4Allocator<RMGOpticalDetectorHit>* RMGOpticalDetectorHitAllocator;

include/RMGOpticalOutputScheme.hh

@@ -18,17 +18,25 @@
 
 #include <vector>
 
+#include "G4AnalysisManager.hh"
+#include "G4GenericMessenger.hh"
+
 #include "RMGVOutputScheme.hh"
 
 class G4Event;
 class RMGOpticalOutputScheme : public RMGVOutputScheme {
 
   public:
 
-    inline RMGOpticalOutputScheme(G4AnalysisManager* ana_man) : RMGVOutputScheme(ana_man) {}
+    RMGOpticalOutputScheme();
 
     void AssignOutputNames(G4AnalysisManager* ana_man) override;
-    void EndOfEventAction(const G4Event*) override;
+    void StoreEvent(const G4Event*) override;
+
+  private:
+
+    std::unique_ptr<G4GenericMessenger> fMessenger;
+    void DefineCommands();
 };
 
 #endif

include/RMGRunAction.hh

@@ -54,7 +54,7 @@ class RMGRunAction : public G4UserRunAction {
       return fOutputDataFields.at(d_type);
     }
     inline void ClearOutputDataFields() {
-      for (auto& el : fOutputDataFields) el.second->clear();
+      for (auto& el : fOutputDataFields) el.second->ClearBeforeEvent();
     }
 
   private:
@@ -66,7 +66,7 @@ class RMGRunAction : public G4UserRunAction {
 
     int fCurrentPrintModulo = -1;
 
-    std::map<RMGHardware::DetectorType, std::unique_ptr<RMGVOutputScheme>> fOutputDataFields;
+    std::map<RMGHardware::DetectorType, std::shared_ptr<RMGVOutputScheme>> fOutputDataFields;
 };
 
 #endif

include/RMGVOutputScheme.hh

@@ -16,7 +16,6 @@
 #ifndef _RMG_V_OUTPUT_SCHEME_HH_
 #define _RMG_V_OUTPUT_SCHEME_HH_
 
-#include <map>
 #include <string>
 
 #include "G4AnalysisManager.hh"
@@ -28,10 +27,16 @@ class RMGVOutputScheme {
 
   public:
 
-    inline RMGVOutputScheme(G4AnalysisManager* ana_man) { this->AssignOutputNames(ana_man); }
-    virtual inline void clear() {};
-    virtual inline void AssignOutputNames(G4AnalysisManager*) {};
-    virtual inline void EndOfEventAction(const G4Event*) {};
+    RMGVOutputScheme() = default;
+
+    // initialization.
+    virtual inline void AssignOutputNames(G4AnalysisManager*) {}
+
+    // functions for individual events.
+    virtual inline void ClearBeforeEvent() {}
+    virtual inline bool ShouldDiscardEvent(const G4Event*) { return false; }
+    virtual inline void StoreEvent(const G4Event*) {}
+
     inline std::string GetNtupleName(int det_uid) { return fmt::format("det{:03}", det_uid); }
 };
 

src/RMGEventAction.cc

@@ -69,8 +69,14 @@ void RMGEventAction::EndOfEventAction(const G4Event* event) {
   auto det_cons = RMGManager::Instance()->GetDetectorConstruction();
   auto active_dets = det_cons->GetActiveDetectorList();
 
+  // Check if any OutputScheme wants to discard this event.
   for (const auto& d_type : active_dets) {
-    fRunAction->GetOutputDataFields(d_type)->EndOfEventAction(event);
+    bool discard_event = fRunAction->GetOutputDataFields(d_type)->ShouldDiscardEvent(event);
+    if (discard_event) return;
+  }
+
+  for (const auto& d_type : active_dets) {
+    fRunAction->GetOutputDataFields(d_type)->StoreEvent(event);
   }
 
   // NOTE: G4analysisManager::AddNtupleRow() must be called here for event-oriented output

src/RMGGermaniumDetector.cc

@@ -42,8 +42,8 @@ G4bool RMGGermaniumDetectorHit::operator==(const RMGGermaniumDetectorHit& right)
 }
 
 void RMGGermaniumDetectorHit::Print() {
-  // TODO: add all fields
   RMGLog::Out(RMGLog::debug, "Detector UID: ", this->detector_uid,
+      " / Particle: ", this->particle_type,
       " / Energy: ", G4BestUnit(this->energy_deposition, "Energy"),
       " / Position: ", this->global_position / CLHEP::m, " m",
       " / Time: ", this->global_time / CLHEP::ns, " ns");
@@ -65,8 +65,6 @@ RMGGermaniumDetector::RMGGermaniumDetector() : G4VSensitiveDetector("Germanium")
   // declare only one hit collection.
   // NOTE: names in the respective output scheme class must match this
   G4VSensitiveDetector::collectionName.insert("Hits");
-
-  this->DefineCommands();
 }
 
 void RMGGermaniumDetector::Initialize(G4HCofThisEvent* hit_coll) {
@@ -134,10 +132,4 @@ bool RMGGermaniumDetector::ProcessHits(G4Step* step, G4TouchableHistory* /*histo
 
 void RMGGermaniumDetector::EndOfEvent(G4HCofThisEvent* /*hit_coll*/) {}
 
-void RMGGermaniumDetector::DefineCommands() {
-
-  fMessenger = std::make_unique<G4GenericMessenger>(this, "/RMG/Detector/Germanium/",
-      "Commands for controlling stuff");
-}
-
 // vim: tabstop=2 shiftwidth=2 expandtab

src/RMGGermaniumOutputScheme.cc

@@ -29,6 +29,8 @@
 
 namespace u = CLHEP;
 
+RMGGermaniumOutputScheme::RMGGermaniumOutputScheme() { this->DefineCommands(); }
+
 // invoked in RMGRunAction::SetupAnalysisManager()
 void RMGGermaniumOutputScheme::AssignOutputNames(G4AnalysisManager* ana_man) {
 
@@ -59,24 +61,58 @@ void RMGGermaniumOutputScheme::AssignOutputNames(G4AnalysisManager* ana_man) {
   }
 }
 
-// invoked in RMGEventAction::EndOfEventAction()
-void RMGGermaniumOutputScheme::EndOfEventAction(const G4Event* event) {
+RMGGermaniumDetectorHitsCollection* RMGGermaniumOutputScheme::GetHitColl(const G4Event* event) {
   auto sd_man = G4SDManager::GetSDMpointer();
 
   auto hit_coll_id = sd_man->GetCollectionID("Germanium/Hits");
   if (hit_coll_id < 0) {
     RMGLog::OutDev(RMGLog::error, "Could not find hit collection Germanium/Hits");
-    return;
+    return nullptr;
   }
 
   auto hit_coll = dynamic_cast<RMGGermaniumDetectorHitsCollection*>(
       event->GetHCofThisEvent()->GetHC(hit_coll_id));
 
   if (!hit_coll) {
     RMGLog::Out(RMGLog::error, "Could not find hit collection associated with event");
-    return;
+    return nullptr;
+  }
+
+  return hit_coll;
+}
+
+// invoked in RMGEventAction::EndOfEventAction()
+bool RMGGermaniumOutputScheme::ShouldDiscardEvent(const G4Event* event) {
+  // exit fast if no threshold is configured.
+  if ((fEdepCutLow < 0 && fEdepCutHigh < 0) || fEdepCutDetectors.empty()) return false;
+
+  auto hit_coll = GetHitColl(event);
+  if (!hit_coll) return false;
+
+  // check defined energy threshold.
+  double event_edep = 0.;
+
+  for (auto hit : *hit_coll->GetVector()) {
+    if (!hit) continue;
+    if (fEdepCutDetectors.find(hit->detector_uid) != fEdepCutDetectors.end())
+      event_edep += hit->energy_deposition;
   }
 
+  if ((fEdepCutLow > 0 && event_edep < fEdepCutLow) ||
+      (fEdepCutHigh > 0 && event_edep > fEdepCutHigh)) {
+    RMGLog::Out(RMGLog::debug, "Discarding event - energy threshold has not been met", event_edep,
+        fEdepCutLow, fEdepCutHigh);
+    return true;
+  }
+
+  return false;
+}
+
+// invoked in RMGEventAction::EndOfEventAction()
+void RMGGermaniumOutputScheme::StoreEvent(const G4Event* event) {
+  auto hit_coll = GetHitColl(event);
+  if (!hit_coll) return;
+
   if (hit_coll->entries() <= 0) {
     RMGLog::OutDev(RMGLog::debug, "Hit collection is empty");
     return;
@@ -109,4 +145,27 @@ void RMGGermaniumOutputScheme::EndOfEventAction(const G4Event* event) {
   }
 }
 
+void RMGGermaniumOutputScheme::DefineCommands() {
+
+  fMessenger = std::make_unique<G4GenericMessenger>(this, "/RMG/Output/Germanium/",
+      "Commands for controlling output from hits in germanium detectors.");
+
+  fMessenger->DeclareMethodWithUnit("SetEdepCutLow", "keV", &RMGGermaniumOutputScheme::SetEdepCutLow)
+      .SetGuidance("Set a lower energy cut that has to be met for this event to be stored.")
+      .SetParameterName("threshold", false)
+      .SetStates(G4State_Idle);
+
+  fMessenger
+      ->DeclareMethodWithUnit("SetEdepCutHigh", "keV", &RMGGermaniumOutputScheme::SetEdepCutHigh)
+      .SetGuidance("Set an upper energy cut that has to be met for this event to be stored.")
+      .SetParameterName("threshold", false)
+      .SetStates(G4State_Idle);
+
+  fMessenger
+      ->DeclareMethod("AddDetectorForEdepThreshold", &RMGGermaniumOutputScheme::AddEdepCutDetector)
+      .SetGuidance("Take this detector into account for the filtering by /EdepThreshold.")
+      .SetParameterName("det_uid", false)
+      .SetStates(G4State_Idle);
+}
+
 // vim: tabstop=2 shiftwidth=2 expandtab

src/RMGHardware.cc

@@ -28,10 +28,12 @@ namespace fs = std::filesystem;
 #include "G4VPhysicalVolume.hh"
 
 #include "RMGGermaniumDetector.hh"
+#include "RMGGermaniumOutputScheme.hh"
 #include "RMGHardwareMessenger.hh"
 #include "RMGLog.hh"
 #include "RMGNavigationTools.hh"
 #include "RMGOpticalDetector.hh"
+#include "RMGOpticalOutputScheme.hh"
 
 #include "magic_enum/magic_enum.hpp"
 
@@ -115,16 +117,24 @@ void RMGHardware::ConstructSDandField() {
           magic_enum::enum_name(v.type));
 
       G4VSensitiveDetector* obj = nullptr;
+      std::shared_ptr<RMGVOutputScheme> output;
       switch (v.type) {
-        case DetectorType::kOptical: obj = new RMGOpticalDetector(); break;
-        case DetectorType::kGermanium: obj = new RMGGermaniumDetector(); break;
+        case DetectorType::kOptical:
+          obj = new RMGOpticalDetector();
+          output = std::make_shared<RMGOpticalOutputScheme>();
+          break;
+        case DetectorType::kGermanium:
+          obj = new RMGGermaniumDetector();
+          output = std::make_shared<RMGGermaniumOutputScheme>();
+          break;
         case DetectorType::kLAr:
         default:
           RMGLog::OutDev(RMGLog::fatal, "No behaviour for sensitive detector type '",
               magic_enum::enum_name<DetectorType>(v.type), "' implemented (implement me)");
       }
       sd_man->AddNewDetector(obj);
       active_dets.emplace(v.type, obj);
+      fActiveOutputSchemes.insert({v.type, output});
     }
 
     // now assign logical volumes to the sensitive detector
@@ -139,14 +149,15 @@ void RMGHardware::ConstructSDandField() {
     RMGLog::OutFormat(RMGLog::debug,
         "Registered new sensitive detector volume of type {}: {} (uid={}, lv={})",
         magic_enum::enum_name(v.type), pv->GetName().c_str(), v.uid, lv->GetName().c_str());
-    fActiveDetectorsInitialized = true;
   }
 
   std::string vec_repr;
   for (const auto& d : fActiveDetectors) vec_repr += std::string(magic_enum::enum_name(d)) + ", ";
   if (vec_repr.size() > 2) vec_repr.erase(vec_repr.size() - 2);
   RMGLog::OutFormat(RMGLog::debug, "List of activated detectors: [{}]", vec_repr);
 
+  fActiveDetectorsInitialized = true;
+
   // copy birks constant from material properties, as it cannot be specified in GDML
   for (G4Material* mat : *G4Material::GetMaterialTable()) {
     auto mpt = mat->GetMaterialPropertiesTable();

src/RMGOpticalDetector.cc

@@ -47,8 +47,6 @@ RMGOpticalDetector::RMGOpticalDetector() : G4VSensitiveDetector("Optical") {
   // declare only one hit collection.
   // NOTE: names in the respective output scheme class must match this
   G4VSensitiveDetector::collectionName.insert("Hits");
-
-  this->DefineCommands();
 }
 
 void RMGOpticalDetector::Initialize(G4HCofThisEvent* hit_coll) {
@@ -117,10 +115,4 @@ bool RMGOpticalDetector::ProcessHits(G4Step* step, G4TouchableHistory* /*history
 
 void RMGOpticalDetector::EndOfEvent(G4HCofThisEvent* /*hit_coll*/) {}
 
-void RMGOpticalDetector::DefineCommands() {
-
-  fMessenger = std::make_unique<G4GenericMessenger>(this, "/RMG/Detector/Optical/",
-      "Commands for controlling stuff");
-}
-
 // vim: tabstop=2 shiftwidth=2 expandtab