add scintillation detector and output scheme

include/RMGHardware.hh

@@ -49,7 +49,7 @@ class RMGHardware : public G4VUserDetectorConstruction {
     enum DetectorType {
       kGermanium,
       kOptical,
-      kLAr
+      kScintillator,
     };
 
     struct DetectorMetadata {

include/RMGScintillatorDetector.hh

@@ -0,0 +1,98 @@
+// 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_SCINTILLATOR_DETECTOR_HH_
+#define _RMG_SCINTILLATOR_DETECTOR_HH_
+
+#include <memory>
+#include <string>
+
+#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<RMGScintillatorDetectorHit>;
+
+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<RMGScintillatorDetectorHit>* RMGScintillatorDetectorHitAllocator;
+
+inline void* RMGScintillatorDetectorHit::operator new(size_t) {
+  if (!RMGScintillatorDetectorHitAllocator)
+    RMGScintillatorDetectorHitAllocator = new G4Allocator<RMGScintillatorDetectorHit>;
+  return (void*)RMGScintillatorDetectorHitAllocator->MallocSingle();
+}
+
+inline void RMGScintillatorDetectorHit::operator delete(void* hit) {
+  RMGScintillatorDetectorHitAllocator->FreeSingle((RMGScintillatorDetectorHit*)hit);
+}
+
+#endif
+
+// vim: tabstop=2 shiftwidth=2 expandtab

include/RMGScintillatorOutputScheme.hh

@@ -0,0 +1,61 @@
+// 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_SCINTILLATOR_OUTPUT_SCHEME_HH_
+#define _RMG_SCINTILLATOR_OUTPUT_SCHEME_HH_
+
+#include <optional>
+#include <set>
+
+#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<G4GenericMessenger> fMessenger;
+    void DefineCommands();
+
+    double fEdepCutLow = -1;
+    double fEdepCutHigh = -1;
+    std::set<int> fEdepCutDetectors;
+};
+
+#endif
+
+// vim: tabstop=2 shiftwidth=2 expandtab

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

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<RMGGermaniumOutputScheme>();
           break;
-        case DetectorType::kLAr:
+        case DetectorType::kScintillator:
+          obj = new RMGScintillatorDetector();
+          output = std::make_shared<RMGScintillatorOutputScheme>();
+          break;
         default:
           RMGLog::OutDev(RMGLog::fatal, "No behaviour for sensitive detector type '",
               magic_enum::enum_name<DetectorType>(v.type), "' implemented (implement me)");

src/RMGScintillatorDetector.cc

@@ -0,0 +1,136 @@
+// 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 "RMGScintillatorDetector.hh"
+
+#include <map>
+#include <stdexcept>
+#include <string>
+
+#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<RMGScintillatorDetectorHit>* 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