From 22736b7878fa3855c8a5932968fba5f4f6188eff Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Fri, 26 Apr 2024 21:50:20 +0200 Subject: [PATCH] add option to discard event with energy threshold in germanium --- include/RMGGermaniumDetector.hh | 3 -- include/RMGGermaniumOutputScheme.hh | 26 +++++++++++- include/RMGHardware.hh | 3 ++ include/RMGOpticalOutputScheme.hh | 6 ++- include/RMGRunAction.hh | 4 +- include/RMGVOutputScheme.hh | 15 ++++--- src/RMGEventAction.cc | 13 +++++- src/RMGGermaniumDetector.cc | 10 +---- src/RMGGermaniumOutputScheme.cc | 64 +++++++++++++++++++++++++++-- src/RMGHardware.cc | 17 ++++++-- src/RMGOpticalOutputScheme.cc | 2 +- src/RMGRunAction.cc | 23 +++-------- 12 files changed, 137 insertions(+), 49 deletions(-) diff --git a/include/RMGGermaniumDetector.hh b/include/RMGGermaniumDetector.hh index 91a8b51..a0211ea 100644 --- a/include/RMGGermaniumDetector.hh +++ b/include/RMGGermaniumDetector.hh @@ -77,9 +77,6 @@ class RMGGermaniumDetector : public G4VSensitiveDetector { private: RMGGermaniumDetectorHitsCollection* fHitsCollection = nullptr; - - std::unique_ptr fMessenger = nullptr; - void DefineCommands(); }; extern G4ThreadLocal G4Allocator* RMGGermaniumDetectorHitAllocator; diff --git a/include/RMGGermaniumOutputScheme.hh b/include/RMGGermaniumOutputScheme.hh index a7c7e2a..1edd818 100644 --- a/include/RMGGermaniumOutputScheme.hh +++ b/include/RMGGermaniumOutputScheme.hh @@ -16,6 +16,12 @@ #ifndef _RMG_GERMANIUM_OUTPUT_SCHEME_HH_ #define _RMG_GERMANIUM_OUTPUT_SCHEME_HH_ +#include + +#include "G4AnalysisManager.hh" +#include "G4GenericMessenger.hh" + +#include "RMGGermaniumDetector.hh" #include "RMGVOutputScheme.hh" class G4Event; @@ -23,10 +29,26 @@ 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 fMessenger; + void DefineCommands(); + + double fEdepCutLow = -1; + double fEdepCutHigh = -1; + std::set fEdepCutDetectors; }; #endif diff --git a/include/RMGHardware.hh b/include/RMGHardware.hh index b1852f4..8a9763e 100644 --- a/include/RMGHardware.hh +++ b/include/RMGHardware.hh @@ -28,6 +28,7 @@ #include "RMGHardwareMessenger.hh" #include "RMGNavigationTools.hh" +#include "RMGVOutputScheme.hh" class G4VPhysicalVolume; class RMGHardware : public G4VUserDetectorConstruction { @@ -62,6 +63,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; } @@ -81,6 +83,7 @@ class RMGHardware : public G4VUserDetectorConstruction { // one element for each sensitive detector physical volume std::map, DetectorMetadata> fDetectorMetadata; std::set fActiveDetectors; + std::map> fActiveOutputSchemes; bool fActiveDetectorsInitialized = false; std::unique_ptr fMessenger; diff --git a/include/RMGOpticalOutputScheme.hh b/include/RMGOpticalOutputScheme.hh index 145f82f..ab7b17f 100644 --- a/include/RMGOpticalOutputScheme.hh +++ b/include/RMGOpticalOutputScheme.hh @@ -18,6 +18,8 @@ #include +#include "G4AnalysisManager.hh" + #include "RMGVOutputScheme.hh" class G4Event; @@ -25,10 +27,10 @@ class RMGOpticalOutputScheme : public RMGVOutputScheme { public: - inline RMGOpticalOutputScheme(G4AnalysisManager* ana_man) : RMGVOutputScheme(ana_man) {} + RMGOpticalOutputScheme() = default; void AssignOutputNames(G4AnalysisManager* ana_man) override; - void EndOfEventAction(const G4Event*) override; + void StoreEvent(const G4Event*) override; }; #endif diff --git a/include/RMGRunAction.hh b/include/RMGRunAction.hh index 622e675..41a31cd 100644 --- a/include/RMGRunAction.hh +++ b/include/RMGRunAction.hh @@ -52,7 +52,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: @@ -62,7 +62,7 @@ class RMGRunAction : public G4UserRunAction { bool fIsAnaManInitialized = false; RMGMasterGenerator* fRMGMasterGenerator = nullptr; - std::map> fOutputDataFields; + std::map> fOutputDataFields; }; #endif diff --git a/include/RMGVOutputScheme.hh b/include/RMGVOutputScheme.hh index 256bd40..3693852 100644 --- a/include/RMGVOutputScheme.hh +++ b/include/RMGVOutputScheme.hh @@ -16,7 +16,6 @@ #ifndef _RMG_V_OUTPUT_SCHEME_HH_ #define _RMG_V_OUTPUT_SCHEME_HH_ -#include #include #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); } }; diff --git a/src/RMGEventAction.cc b/src/RMGEventAction.cc index 0fa76f9..86e1cfe 100644 --- a/src/RMGEventAction.cc +++ b/src/RMGEventAction.cc @@ -69,8 +69,19 @@ 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. + bool discard_event = false; for (const auto& d_type : active_dets) { - fRunAction->GetOutputDataFields(d_type)->EndOfEventAction(event); + discard_event &= fRunAction->GetOutputDataFields(d_type)->ShouldDiscardEvent(event); + } + + if (discard_event) { + RMGLog::Out(RMGLog::debug, "Discarding event as requested by output scheme."); + 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 diff --git a/src/RMGGermaniumDetector.cc b/src/RMGGermaniumDetector.cc index 1cdc6e6..569ab90 100644 --- a/src/RMGGermaniumDetector.cc +++ b/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(this, "/RMG/Detector/Germanium/", - "Commands for controlling stuff"); -} - // vim: tabstop=2 shiftwidth=2 expandtab diff --git a/src/RMGGermaniumOutputScheme.cc b/src/RMGGermaniumOutputScheme.cc index 1348d16..be449ac 100644 --- a/src/RMGGermaniumOutputScheme.cc +++ b/src/RMGGermaniumOutputScheme.cc @@ -27,6 +27,8 @@ namespace u = CLHEP; +RMGGermaniumOutputScheme::RMGGermaniumOutputScheme() { this->DefineCommands(); } + // invoked in RMGRunAction::SetupAnalysisManager() void RMGGermaniumOutputScheme::AssignOutputNames(G4AnalysisManager* ana_man) { @@ -52,14 +54,13 @@ 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( @@ -67,9 +68,41 @@ void RMGGermaniumOutputScheme::EndOfEventAction(const G4Event* event) { 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) { + if ((fEdepCutLow < 0 && fEdepCutHigh) || !fEdepCutDetectors.empty()) return false; + + auto hit_coll = GetHitColl(event); + if (!hit_coll) return false; + + // check defined energy threshold. + auto event_edep = 0; + + for (auto hit : *hit_coll->GetVector()) { + if (!hit) continue; + auto search = fEdepCutDetectors.find(hit->detector_uid); + if (search != fEdepCutDetectors.end()) event_edep += hit->energy_deposition; } + if (event_edep < fEdepCutLow || event_edep > fEdepCutHigh) { + RMGLog::Out(RMGLog::debug, "Discarding event - energy threshold has not been met"); + 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; @@ -102,4 +135,27 @@ void RMGGermaniumOutputScheme::EndOfEventAction(const G4Event* event) { } } +void RMGGermaniumOutputScheme::DefineCommands() { + + fMessenger = std::make_unique(this, "/RMG/Output/Germanium/", + "Commands for controlling output from hits in germanium detectors."); + + fMessenger->DeclareMethodWithUnit("SetEdepCutLow", "keV", &RMGGermaniumOutputScheme::SetEdepCutLow) + .SetGuidance("Set an energy threshold 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 energy threshold 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 diff --git a/src/RMGHardware.cc b/src/RMGHardware.cc index b86418a..4e66a39 100644 --- a/src/RMGHardware.cc +++ b/src/RMGHardware.cc @@ -27,10 +27,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" @@ -108,9 +110,16 @@ void RMGHardware::ConstructSDandField() { magic_enum::enum_name(v.type)); G4VSensitiveDetector* obj = nullptr; + std::shared_ptr 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(); + break; + case DetectorType::kGermanium: + obj = new RMGGermaniumDetector(); + output = std::make_shared(); + break; case DetectorType::kLAr: default: RMGLog::OutDev(RMGLog::fatal, "No behaviour for sensitive detector type '", @@ -118,6 +127,7 @@ void RMGHardware::ConstructSDandField() { } sd_man->AddNewDetector(obj); active_dets.emplace(v.type, obj); + fActiveOutputSchemes.insert({v.type, output}); } // now assign logical volumes to the sensitive detector @@ -132,7 +142,6 @@ 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; @@ -140,6 +149,8 @@ void RMGHardware::ConstructSDandField() { 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(); diff --git a/src/RMGOpticalOutputScheme.cc b/src/RMGOpticalOutputScheme.cc index 192c1a6..e8e47ce 100644 --- a/src/RMGOpticalOutputScheme.cc +++ b/src/RMGOpticalOutputScheme.cc @@ -49,7 +49,7 @@ void RMGOpticalOutputScheme::AssignOutputNames(G4AnalysisManager* ana_man) { } // invoked in RMGEventAction::EndOfEventAction() -void RMGOpticalOutputScheme::EndOfEventAction(const G4Event* event) { +void RMGOpticalOutputScheme::StoreEvent(const G4Event* event) { auto sd_man = G4SDManager::GetSDMpointer(); auto hit_coll_id = sd_man->GetCollectionID("Optical/Hits"); diff --git a/src/RMGRunAction.cc b/src/RMGRunAction.cc index 938760f..4760416 100644 --- a/src/RMGRunAction.cc +++ b/src/RMGRunAction.cc @@ -71,24 +71,13 @@ void RMGRunAction::SetupAnalysisManager() { auto det_cons = RMGManager::Instance()->GetDetectorConstruction(); for (const auto& d_type : det_cons->GetActiveDetectorList()) { - RMGLog::OutFormatDev(RMGLog::debug, - "Initializing output scheme for sensitive detector type '{}'", magic_enum::enum_name(d_type)); - - // instantiate concrete output scheme class - // TODO: allow for library user to register own output scheme - switch (d_type) { - case RMGHardware::kOptical: - fOutputDataFields.emplace(d_type, std::make_unique(ana_man)); - fOutputDataFields[d_type]->AssignOutputNames(ana_man); - break; - case RMGHardware::kGermanium: - fOutputDataFields.emplace(d_type, std::make_unique(ana_man)); - fOutputDataFields[d_type]->AssignOutputNames(ana_man); - break; - default: - RMGLog::OutDev(RMGLog::fatal, "No output scheme sensitive detector type '", - magic_enum::enum_name(d_type), "' implemented (implement me)"); + fOutputDataFields[d_type] = det_cons->GetActiveOutputScheme(d_type); + if (!fOutputDataFields[d_type]) { + RMGLog::OutDev(RMGLog::fatal, "No output scheme sensitive detector type '", + magic_enum::enum_name(d_type), "' implemented (implement me)"); } + + fOutputDataFields[d_type]->AssignOutputNames(ana_man); } }