From 71d109c4d7c7447f87dc91a7255b5a512ae909a6 Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Fri, 19 Jul 2024 23:37:47 +0200 Subject: [PATCH] add scintillation detector and output scheme --- include/RMGHardware.hh | 2 +- include/RMGScintillatorDetector.hh | 98 +++++++++++++ include/RMGScintillatorOutputScheme.hh | 61 ++++++++ src/CMakeLists.txt | 4 + src/RMGHardware.cc | 7 +- src/RMGScintillatorDetector.cc | 136 +++++++++++++++++ src/RMGScintillatorOutputScheme.cc | 193 +++++++++++++++++++++++++ 7 files changed, 499 insertions(+), 2 deletions(-) create mode 100644 include/RMGScintillatorDetector.hh create mode 100644 include/RMGScintillatorOutputScheme.hh create mode 100644 src/RMGScintillatorDetector.cc create mode 100644 src/RMGScintillatorOutputScheme.cc diff --git a/include/RMGHardware.hh b/include/RMGHardware.hh index cfe7b71..9fe2367 100644 --- a/include/RMGHardware.hh +++ b/include/RMGHardware.hh @@ -49,7 +49,7 @@ class RMGHardware : public G4VUserDetectorConstruction { enum DetectorType { kGermanium, kOptical, - kLAr + kScintillator, }; struct DetectorMetadata { diff --git a/include/RMGScintillatorDetector.hh b/include/RMGScintillatorDetector.hh new file mode 100644 index 0000000..0fa4576 --- /dev/null +++ b/include/RMGScintillatorDetector.hh @@ -0,0 +1,98 @@ +// Copyright (C) 2022 Luigi Pertoldi +// +// 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 . + +#ifndef _RMG_SCINTILLATOR_DETECTOR_HH_ +#define _RMG_SCINTILLATOR_DETECTOR_HH_ + +#include +#include + +#include "G4Allocator.hh" +#include "G4THitsCollection.hh" +#include "G4ThreeVector.hh" +#include "G4VHit.hh" +#include "G4VSensitiveDetector.hh" + +class RMGScintillatorDetectorHit : public G4VHit { + + public: + + RMGScintillatorDetectorHit() = default; + ~RMGScintillatorDetectorHit() = default; + + RMGScintillatorDetectorHit(RMGScintillatorDetectorHit const&) = delete; + RMGScintillatorDetectorHit& operator=(RMGScintillatorDetectorHit const&) = delete; + RMGScintillatorDetectorHit(RMGScintillatorDetectorHit&&) = delete; + RMGScintillatorDetectorHit& operator=(RMGScintillatorDetectorHit&&) = delete; + + bool operator==(const RMGScintillatorDetectorHit&) const; + + inline void* operator new(size_t); + inline void operator delete(void*); + + void Print() override; + void Draw() override; + + int detector_uid = -1; + int particle_type = -1; + float energy_deposition = -1; + G4ThreeVector global_position_pre; + G4ThreeVector global_position_post; + double global_time = -1; + double velocity_pre = -1; + double velocity_post = -1; +}; + +using RMGScintillatorDetectorHitsCollection = G4THitsCollection; + +class G4Step; +class G4HCofThisEvent; +class G4TouchableHistory; +class RMGScintillatorDetector : public G4VSensitiveDetector { + + public: + + RMGScintillatorDetector(); + ~RMGScintillatorDetector() = default; + + RMGScintillatorDetector(RMGScintillatorDetector const&) = delete; + RMGScintillatorDetector& operator=(RMGScintillatorDetector const&) = delete; + RMGScintillatorDetector(RMGScintillatorDetector&&) = delete; + RMGScintillatorDetector& operator=(RMGScintillatorDetector&&) = delete; + + void Initialize(G4HCofThisEvent* hit_coll) override; + bool ProcessHits(G4Step* step, G4TouchableHistory* history) override; + void EndOfEvent(G4HCofThisEvent* hit_coll) override; + + private: + + RMGScintillatorDetectorHitsCollection* fHitsCollection = nullptr; +}; + +extern G4ThreadLocal G4Allocator* RMGScintillatorDetectorHitAllocator; + +inline void* RMGScintillatorDetectorHit::operator new(size_t) { + if (!RMGScintillatorDetectorHitAllocator) + RMGScintillatorDetectorHitAllocator = new G4Allocator; + return (void*)RMGScintillatorDetectorHitAllocator->MallocSingle(); +} + +inline void RMGScintillatorDetectorHit::operator delete(void* hit) { + RMGScintillatorDetectorHitAllocator->FreeSingle((RMGScintillatorDetectorHit*)hit); +} + +#endif + +// vim: tabstop=2 shiftwidth=2 expandtab diff --git a/include/RMGScintillatorOutputScheme.hh b/include/RMGScintillatorOutputScheme.hh new file mode 100644 index 0000000..dc9ca99 --- /dev/null +++ b/include/RMGScintillatorOutputScheme.hh @@ -0,0 +1,61 @@ +// Copyright (C) 2022 Luigi Pertoldi +// +// 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 . + +#ifndef _RMG_SCINTILLATOR_OUTPUT_SCHEME_HH_ +#define _RMG_SCINTILLATOR_OUTPUT_SCHEME_HH_ + +#include +#include + +#include "G4AnalysisManager.hh" +#include "G4GenericMessenger.hh" + +#include "RMGScintillatorDetector.hh" +#include "RMGVOutputScheme.hh" + +class G4Event; +class RMGScintillatorOutputScheme : public RMGVOutputScheme { + + public: + + RMGScintillatorOutputScheme(); + + void AssignOutputNames(G4AnalysisManager* ana_man) 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); } + + protected: + + [[nodiscard]] inline std::string GetNtuplenameFlat() const override { return "scintillator"; } + + private: + + RMGScintillatorDetectorHitsCollection* GetHitColl(const G4Event*); + + std::unique_ptr fMessenger; + void DefineCommands(); + + double fEdepCutLow = -1; + double fEdepCutHigh = -1; + std::set fEdepCutDetectors; +}; + +#endif + +// vim: tabstop=2 shiftwidth=2 expandtab diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ab360b9..079443f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -22,6 +22,8 @@ set(PROJECT_PUBLIC_HEADERS ${_root}/include/RMGPhysics.hh ${_root}/include/RMGRunAction.hh ${_root}/include/RMGRun.hh + ${_root}/include/RMGScintillatorDetector.hh + ${_root}/include/RMGScintillatorOutputScheme.hh ${_root}/include/RMGStackingAction.hh ${_root}/include/RMGSteppingAction.hh ${_root}/include/RMGTools.hh @@ -52,6 +54,8 @@ set(PROJECT_SOURCES ${_root}/src/RMGOpticalOutputScheme.cc ${_root}/src/RMGPhysics.cc ${_root}/src/RMGRunAction.cc + ${_root}/src/RMGScintillatorDetector.cc + ${_root}/src/RMGScintillatorOutputScheme.cc ${_root}/src/RMGStackingAction.cc ${_root}/src/RMGSteppingAction.cc ${_root}/src/RMGTrackingAction.cc diff --git a/src/RMGHardware.cc b/src/RMGHardware.cc index ee8e1fc..0c45dbd 100644 --- a/src/RMGHardware.cc +++ b/src/RMGHardware.cc @@ -34,6 +34,8 @@ namespace fs = std::filesystem; #include "RMGNavigationTools.hh" #include "RMGOpticalDetector.hh" #include "RMGOpticalOutputScheme.hh" +#include "RMGScintillatorDetector.hh" +#include "RMGScintillatorOutputScheme.hh" #include "RMGVertexOutputScheme.hh" #include "magic_enum/magic_enum.hpp" @@ -133,7 +135,10 @@ void RMGHardware::ConstructSDandField() { obj = new RMGGermaniumDetector(); output = std::make_shared(); break; - case DetectorType::kLAr: + case DetectorType::kScintillator: + obj = new RMGScintillatorDetector(); + output = std::make_shared(); + break; default: RMGLog::OutDev(RMGLog::fatal, "No behaviour for sensitive detector type '", magic_enum::enum_name(v.type), "' implemented (implement me)"); diff --git a/src/RMGScintillatorDetector.cc b/src/RMGScintillatorDetector.cc new file mode 100644 index 0000000..df240cd --- /dev/null +++ b/src/RMGScintillatorDetector.cc @@ -0,0 +1,136 @@ +// Copyright (C) 2022 Luigi Pertoldi +// +// 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 . + +#include "RMGScintillatorDetector.hh" + +#include +#include +#include + +#include "G4Circle.hh" +#include "G4HCofThisEvent.hh" +#include "G4OpticalPhoton.hh" +#include "G4SDManager.hh" +#include "G4Step.hh" +#include "G4UnitsTable.hh" +#include "G4VVisManager.hh" + +#include "RMGHardware.hh" +#include "RMGLog.hh" +#include "RMGManager.hh" + +G4ThreadLocal G4Allocator* RMGScintillatorDetectorHitAllocator = nullptr; + +// NOTE: does this make sense? +G4bool RMGScintillatorDetectorHit::operator==(const RMGScintillatorDetectorHit& right) const { + return this == &right; +} + +void RMGScintillatorDetectorHit::Print() { + RMGLog::Out(RMGLog::debug, "Detector UID: ", this->detector_uid, + " / Particle: ", this->particle_type, + " / Energy: ", G4BestUnit(this->energy_deposition, "Energy"), + " / Position: ", this->global_position_pre / CLHEP::m, " m", + " / Time: ", this->global_time / CLHEP::ns, " ns"); +} + +void RMGScintillatorDetectorHit::Draw() { + const auto vis_man = G4VVisManager::GetConcreteInstance(); + if (vis_man and this->energy_deposition > 0) { + G4Circle circle(this->global_position_pre); + circle.SetScreenSize(5); + circle.SetFillStyle(G4Circle::filled); + circle.SetVisAttributes(G4VisAttributes(G4Colour(1, 0, 0))); + vis_man->Draw(circle); + } +} + +RMGScintillatorDetector::RMGScintillatorDetector() : G4VSensitiveDetector("Scintillator") { + + // declare only one hit collection. + // NOTE: names in the respective output scheme class must match this + G4VSensitiveDetector::collectionName.insert("Hits"); +} + +void RMGScintillatorDetector::Initialize(G4HCofThisEvent* hit_coll) { + + // create hits collection object + // NOTE: assumes there is only one collection name (see constructor) + fHitsCollection = + new RMGScintillatorDetectorHitsCollection(G4VSensitiveDetector::SensitiveDetectorName, + G4VSensitiveDetector::collectionName[0]); + + // associate it with the G4HCofThisEvent object + auto hc_id = G4SDManager::GetSDMpointer()->GetCollectionID(G4VSensitiveDetector::collectionName[0]); + hit_coll->AddHitsCollection(hc_id, fHitsCollection); +} + +bool RMGScintillatorDetector::ProcessHits(G4Step* step, G4TouchableHistory* /*history*/) { + + RMGLog::OutDev(RMGLog::debug, "Processing scintillator detector hits"); + + // return if no energy is deposited + if (step->GetTotalEnergyDeposit() == 0) return false; + // ignore optical photons + if (step->GetTrack()->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) return false; + + // we're going to use info from the pre-step point + const auto prestep = step->GetPreStepPoint(); + const auto poststep = step->GetPostStepPoint(); + + // locate us + const auto pv_name = prestep->GetTouchableHandle()->GetVolume()->GetName(); + const auto pv_copynr = prestep->GetTouchableHandle()->GetCopyNumber(); + + // check if physical volume is registered as sermanium detector + const auto det_cons = RMGManager::Instance()->GetDetectorConstruction(); + try { + auto d_type = det_cons->GetDetectorMetadata({pv_name, pv_copynr}).type; + if (d_type != RMGHardware::kScintillator) { + RMGLog::OutFormatDev(RMGLog::debug, + "Volume '{}' (copy nr. {} not registered as scintillator detector", pv_name, pv_copynr); + return false; + } + } catch (const std::out_of_range& e) { + RMGLog::OutFormatDev(RMGLog::debug, "Volume '{}' (copy nr. {}) not registered as detector", + pv_name, pv_copynr); + return false; + } + + // retrieve unique id for persistency + auto det_uid = det_cons->GetDetectorMetadata({pv_name, pv_copynr}).uid; + + RMGLog::OutDev(RMGLog::debug, "Hit in scintillator detector nr. ", det_uid, " detected"); + + // create a new hit and fill it + auto* hit = new RMGScintillatorDetectorHit(); + hit->detector_uid = det_uid; + hit->particle_type = step->GetTrack()->GetDefinition()->GetPDGEncoding(); + hit->energy_deposition = step->GetTotalEnergyDeposit(); + hit->global_position_pre = prestep->GetPosition(); + hit->global_position_post = poststep->GetPosition(); + hit->global_time = prestep->GetGlobalTime(); + hit->velocity_pre = prestep->GetVelocity(); + hit->velocity_post = poststep->GetVelocity(); + + // register the hit in the hit collection for the event + fHitsCollection->insert(hit); + + return true; +} + +void RMGScintillatorDetector::EndOfEvent(G4HCofThisEvent* /*hit_coll*/) {} + +// vim: tabstop=2 shiftwidth=2 expandtab diff --git a/src/RMGScintillatorOutputScheme.cc b/src/RMGScintillatorOutputScheme.cc new file mode 100644 index 0000000..fb6ab29 --- /dev/null +++ b/src/RMGScintillatorOutputScheme.cc @@ -0,0 +1,193 @@ +// Copyright (C) 2022 Luigi Pertoldi +// +// 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 . + +#include "RMGScintillatorOutputScheme.hh" + +#include + +#include "G4AnalysisManager.hh" +#include "G4Event.hh" +#include "G4HCtable.hh" +#include "G4OpticalPhoton.hh" +#include "G4SDManager.hh" + +#include "RMGHardware.hh" +#include "RMGLog.hh" +#include "RMGManager.hh" +#include "RMGScintillatorDetector.hh" + +namespace u = CLHEP; + +RMGScintillatorOutputScheme::RMGScintillatorOutputScheme() { this->DefineCommands(); } + +// invoked in RMGRunAction::SetupAnalysisManager() +void RMGScintillatorOutputScheme::AssignOutputNames(G4AnalysisManager* ana_man) { + + auto rmg_man = RMGManager::Instance(); + const auto det_cons = rmg_man->GetDetectorConstruction(); + const auto detectors = det_cons->GetDetectorMetadataMap(); + + std::set registered_uids; + std::set registered_ntuples; + for (auto&& det : detectors) { + if (det.second.type != RMGHardware::kScintillator) continue; + + // do not register the ntuple twice if two detectors share their uid. + auto had_uid = registered_uids.emplace(det.second.uid); + if (!had_uid.second) continue; + + auto ntuple_name = this->GetNtupleName(det.second.uid); + auto had_name = registered_ntuples.emplace(ntuple_name); + if (!had_name.second) continue; + + auto id = + rmg_man->RegisterNtuple(det.second.uid, ana_man->CreateNtuple(ntuple_name, "Event data")); + + ana_man->CreateNtupleIColumn(id, "evtid"); + if (!fNtuplePerDetector) { ana_man->CreateNtupleIColumn(id, "det_uid"); } + ana_man->CreateNtupleIColumn(id, "particle"); + ana_man->CreateNtupleDColumn(id, "edep"); + ana_man->CreateNtupleDColumn(id, "time"); + ana_man->CreateNtupleDColumn(id, "xloc_pre"); + ana_man->CreateNtupleDColumn(id, "yloc_pre"); + ana_man->CreateNtupleDColumn(id, "zloc_pre"); + ana_man->CreateNtupleDColumn(id, "xloc_post"); + ana_man->CreateNtupleDColumn(id, "yloc_post"); + ana_man->CreateNtupleDColumn(id, "zloc_post"); + ana_man->CreateNtupleDColumn(id, "v_pre"); + ana_man->CreateNtupleDColumn(id, "v_post"); + + ana_man->FinishNtuple(id); + } +} + +RMGScintillatorDetectorHitsCollection* RMGScintillatorOutputScheme::GetHitColl(const G4Event* event) { + auto sd_man = G4SDManager::GetSDMpointer(); + + auto hit_coll_id = sd_man->GetCollectionID("Scintillator/Hits"); + if (hit_coll_id < 0) { + RMGLog::OutDev(RMGLog::error, "Could not find hit collection Scintillator/Hits"); + return nullptr; + } + + auto hit_coll = dynamic_cast( + event->GetHCofThisEvent()->GetHC(hit_coll_id)); + + if (!hit_coll) { + RMGLog::Out(RMGLog::error, "Could not find hit collection associated with event"); + return nullptr; + } + + return hit_coll; +} + +// invoked in RMGEventAction::EndOfEventAction() +bool RMGScintillatorOutputScheme::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 RMGScintillatorOutputScheme::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; + } else { + RMGLog::OutDev(RMGLog::debug, "Hit collection contains ", hit_coll->entries(), " hits"); + } + + auto rmg_man = RMGManager::Instance(); + if (rmg_man->IsPersistencyEnabled()) { + RMGLog::OutDev(RMGLog::debug, "Filling persistent data vectors"); + const auto ana_man = G4AnalysisManager::Instance(); + + for (auto hit : *hit_coll->GetVector()) { + if (!hit) continue; + hit->Print(); + + auto ntupleid = rmg_man->GetNtupleID(hit->detector_uid); + + int col_id = 0; + ana_man->FillNtupleIColumn(ntupleid, col_id++, event->GetEventID()); + if (!fNtuplePerDetector) { + ana_man->FillNtupleIColumn(ntupleid, col_id++, hit->detector_uid); + } + ana_man->FillNtupleIColumn(ntupleid, col_id++, hit->particle_type); + ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->energy_deposition / u::keV); + ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->global_time / u::ns); + ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->global_position_pre.getX() / u::m); + ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->global_position_pre.getY() / u::m); + ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->global_position_pre.getZ() / u::m); + ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->global_position_post.getX() / u::m); + ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->global_position_post.getY() / u::m); + ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->global_position_post.getZ() / u::m); + ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->velocity_pre); + ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->velocity_post); + + // NOTE: must be called here for hit-oriented output + ana_man->AddNtupleRow(ntupleid); + } + } +} + +void RMGScintillatorOutputScheme::DefineCommands() { + + fMessenger = std::make_unique(this, "/RMG/Output/Scintillator/", + "Commands for controlling output from hits in scintillator detectors."); + + fMessenger + ->DeclareMethodWithUnit("SetEdepCutLow", "keV", &RMGScintillatorOutputScheme::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", &RMGScintillatorOutputScheme::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", &RMGScintillatorOutputScheme::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