Add optical SD check and output in post boundary step

src/celeritas/ext/GeantSetup.cc

@@ -11,6 +11,7 @@
 #include <G4ParticleTable.hh>
 #include <G4RunManager.hh>
 #include <G4VPhysicalVolume.hh>
+#include <G4VSensitiveDetector.hh>
 #include <G4VUserDetectorConstruction.hh>
 #include <G4Version.hh>
 #if G4VERSION_NUMBER >= 1070
@@ -30,29 +31,26 @@
 #include "geocel/GeantUtils.hh"
 #include "geocel/ScopedGeantExceptionHandler.hh"
 #include "geocel/ScopedGeantLogger.hh"
+#include "celeritas/g4/DetectorConstruction.hh"
 
 #include "EmPhysicsList.hh"
 
 namespace celeritas
 {
 namespace
 {
+
 //---------------------------------------------------------------------------//
-/*!
- * Load the detector geometry from a GDML input file.
- */
-class DetectorConstruction : public G4VUserDetectorConstruction
+//! Placeholder SD class for generating model data from GDML
+class GdmlSensitiveDetector final : public G4VSensitiveDetector
 {
   public:
-    explicit DetectorConstruction(G4VPhysicalVolume* world) : world_{world}
+    GdmlSensitiveDetector(std::string const& name) : G4VSensitiveDetector{name}
     {
-        CELER_ENSURE(world_);
     }
 
-    G4VPhysicalVolume* Construct() override { return world_; }
-
-  private:
-    G4VPhysicalVolume* world_;
+    void Initialize(G4HCofThisEvent*) final {}
+    bool ProcessHits(G4Step*, G4TouchableHistory*) final { return false; }
 };
 
 //---------------------------------------------------------------------------//
@@ -102,17 +100,21 @@ GeantSetup::GeantSetup(std::string const& gdml_filename, Options options)
     ScopedGeantLogger scoped_logger(celeritas::world_logger());
     ScopedGeantExceptionHandler scoped_exceptions;
 
-    G4VPhysicalVolume* world{nullptr};
+    DetectorConstruction const* p_detector = nullptr;
     {
         CELER_LOG(status) << "Initializing Geant4 geometry and physics list";
 
         // Load GDML and reference the world pointer
         // TODO: pass GdmlLoader options through SetupOptions
-        world = load_gdml(gdml_filename);
-        CELER_ASSERT(world);
 
         // Construct the geometry
-        auto detector = std::make_unique<DetectorConstruction>(world);
+        // auto detector = std::make_unique<DetectorConstruction>(world);
+        auto build_sd_functor = [](std::string const& name) {
+            return std::make_unique<GdmlSensitiveDetector>(name);
+        };
+        auto detector = std::make_unique<DetectorConstruction>(
+            gdml_filename, build_sd_functor);
+        p_detector = detector.get();
         run_manager_->SetUserInitialization(detector.release());
 
         // Construct the physics
@@ -131,7 +133,8 @@ GeantSetup::GeantSetup(std::string const& gdml_filename, Options options)
 
     {
         // Create non-owning Geant4 geo wrapper and save as tracking geometry
-        geo_ = std::make_shared<GeantGeoParams>(world, Ownership::reference);
+        geo_ = std::make_shared<GeantGeoParams>(p_detector->world(),
+                                                Ownership::reference);
         celeritas::global_geant_geo(geo_);
     }
 

src/celeritas/global/CoreParams.cc

@@ -54,6 +54,7 @@
 #include "celeritas/track/SimParams.hh"  // IWYU pragma: keep
 #include "celeritas/track/SortTracksAction.hh"
 #include "celeritas/track/TrackInitParams.hh"  // IWYU pragma: keep
+#include "celeritas/user/SDParams.hh"
 
 #include "ActionInterface.hh"
 
@@ -91,6 +92,7 @@ build_params_refs(CoreParams::Input const& p, CoreScalars const& scalars)
     ref.cutoffs = get_ref<M>(*p.cutoff);
     ref.physics = get_ref<M>(*p.physics);
     ref.rng = get_ref<M>(*p.rng);
+    ref.detector = get_ref<M>(*p.detector);
     ref.sim = get_ref<M>(*p.sim);
     ref.surface = get_ref<M>(*p.surface);
     ref.init = get_ref<M>(*p.init);

src/celeritas/global/CoreParams.hh

@@ -28,6 +28,7 @@ class MpiCommunicator;
 class OutputRegistry;
 class ParticleParams;
 class PhysicsParams;
+class SDParams;
 class SimParams;
 class SurfaceParams;
 class TrackInitParams;
@@ -56,6 +57,7 @@ class CoreParams final : public ParamsDataInterface<CoreParamsData>
     using SPConstCutoff = std::shared_ptr<CutoffParams const>;
     using SPConstPhysics = std::shared_ptr<PhysicsParams const>;
     using SPConstRng = std::shared_ptr<RngParams const>;
+    using SPConstSensDet = std::shared_ptr<SDParams const>;
     using SPConstSim = std::shared_ptr<SimParams const>;
     using SPConstSurface = std::shared_ptr<SurfaceParams const>;
     using SPConstTrackInit = std::shared_ptr<TrackInitParams const>;
@@ -83,6 +85,7 @@ class CoreParams final : public ParamsDataInterface<CoreParamsData>
         SPConstCutoff cutoff;
         SPConstPhysics physics;
         SPConstRng rng;
+        SPConstSensDet detector;
         SPConstSim sim;
         SPConstSurface surface;
         SPConstTrackInit init;
@@ -125,6 +128,7 @@ class CoreParams final : public ParamsDataInterface<CoreParamsData>
     SPConstCutoff const& cutoff() const { return input_.cutoff; }
     SPConstPhysics const& physics() const { return input_.physics; }
     SPConstRng const& rng() const { return input_.rng; }
+    SPConstSensDet const& detector() const { return input_.detector; }
     SPConstSim const& sim() const { return input_.sim; }
     SPConstSurface const& surface() const { return input_.surface; }
     SPConstTrackInit const& init() const { return input_.init; }

src/celeritas/global/CoreTrackData.hh

@@ -21,6 +21,7 @@
 #include "celeritas/phys/PhysicsData.hh"
 #include "celeritas/track/SimData.hh"
 #include "celeritas/track/TrackInitData.hh"
+#include "celeritas/user/SDData.hh"
 
 #include "CoreTrackDataFwd.hh"
 
@@ -72,6 +73,7 @@ struct CoreParamsData
     CutoffParamsData<W, M> cutoffs;
     PhysicsParamsData<W, M> physics;
     RngParamsData<W, M> rng;
+    SDParamsData<W, M> detector;
     SimParamsData<W, M> sim;
     SurfaceParamsData<W, M> surface;
     TrackInitParamsData<W, M> init;
@@ -98,6 +100,7 @@ struct CoreParamsData
         cutoffs = other.cutoffs;
         physics = other.physics;
         rng = other.rng;
+        detector = other.detector;
         sim = other.sim;
         surface = other.surface;
         init = other.init;

src/celeritas/optical/CoreParams.cc

@@ -52,7 +52,7 @@ build_params_refs(CoreParams::Input const& p, CoreScalars const& scalars)
     ref.surface = get_ref<M>(*p.surface);
     ref.surface_physics = get_ref<M>(*p.surface_physics);
     ref.rng = get_ref<M>(*p.rng);
-    // TODO: Get detectors ref
+    ref.detectors = get_ref<M>(*p.detectors);
     if (p.cherenkov)
     {
         ref.cherenkov = get_ref<M>(*p.cherenkov);
@@ -119,17 +119,11 @@ CoreParams::CoreParams(Input&& input) : input_(std::move(input))
     CP_VALIDATE_INPUT(action_reg);
     CP_VALIDATE_INPUT(gen_reg);
     CP_VALIDATE_INPUT(max_streams);
+    CP_VALIDATE_INPUT(detectors);
 #undef CP_VALIDATE_INPUT
 
     CELER_EXPECT(input_);
 
-    // TODO: provide detectors in input, passing from core params
-    detectors_ = input_.detectors;
-    if (!detectors_)
-    {
-        detectors_ = std::make_shared<SDParams>();
-    }
-
     ScopedMem record_mem("optical::CoreParams.construct");
 
     // Construct always-on actions and save their IDs

src/celeritas/optical/CoreParams.hh

@@ -112,7 +112,7 @@ class CoreParams final : public ParamsDataInterface<CoreParamsData>
     }
     SPActionRegistry const& action_reg() const { return input_.action_reg; }
     SPGeneratorRegistry const& gen_reg() const { return input_.gen_reg; }
-    SPConstDetectors const& detectors() const { return detectors_; }
+    SPConstDetectors const& detectors() const { return input_.detectors; }
     SPConstCherenkov const& cherenkov() const { return input_.cherenkov; }
     SPConstScintillation const& scintillation() const
     {
@@ -137,8 +137,6 @@ class CoreParams final : public ParamsDataInterface<CoreParamsData>
 
     // Copy of DeviceRef in device memory
     DeviceVector<DeviceRef> device_ref_vec_;
-
-    SPConstDetectors detectors_;
 };
 
 //---------------------------------------------------------------------------//

src/celeritas/optical/CoreTrackData.hh

@@ -12,6 +12,7 @@
 #include "geocel/SurfaceData.hh"
 #include "celeritas/Types.hh"
 #include "celeritas/geo/GeoData.hh"
+#include "celeritas/user/SDData.hh"
 
 #include "CoreTrackDataFwd.hh"
 #include "MaterialData.hh"
@@ -55,6 +56,7 @@ struct CoreParamsData
     RngParamsData<W, M> rng;
     SurfaceParamsData<W, M> surface;
     SurfacePhysicsParamsData<W, M> surface_physics;
+    SDParamsData<W, M> detectors;
     CherenkovData<W, M> cherenkov;
     ScintillationData<W, M> scintillation;
 

src/celeritas/optical/CoreTrackView.hh

@@ -14,6 +14,7 @@
 #include "MaterialView.hh"
 #include "ParticleTrackView.hh"
 #include "PhysicsTrackView.hh"
+#include "SensitiveDetectorView.hh"
 #include "SimTrackView.hh"
 #include "TrackInitializer.hh"
 #include "surface/SurfacePhysicsTrackView.hh"
@@ -78,6 +79,9 @@ class CoreTrackView
     // Return a surface physics view
     inline CELER_FUNCTION SurfacePhysicsTrackView surface_physics() const;
 
+    // Return a sensitive detector view
+    inline CELER_FUNCTION SensitiveDetectorView sensitive_detectors() const;
+
     // Return an RNG engine
     inline CELER_FUNCTION RngEngine rng() const;
 
@@ -246,6 +250,16 @@ CELER_FUNCTION auto CoreTrackView::surface_physics() const
                                    this->track_slot_id()};
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Return a sensitive detector view.
+ */
+CELER_FUNCTION auto CoreTrackView::sensitive_detectors() const
+    -> SensitiveDetectorView
+{
+    return SensitiveDetectorView{params_.detectors};
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Return the RNG engine.

src/celeritas/optical/SensitiveDetectorView.hh

@@ -0,0 +1,44 @@
+//------------------------------- -*- C++ -*- -------------------------------//
+// Copyright Celeritas contributors: see top-level COPYRIGHT file for details
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/optical/SensitiveDetectorView.hh
+//---------------------------------------------------------------------------//
+#pragma once
+
+#include "geocel/Types.hh"
+#include "celeritas/Quantities.hh"
+#include "celeritas/user/SDData.hh"
+
+namespace celeritas
+{
+namespace optical
+{
+
+class SensitiveDetectorView
+{
+    using SDParamsRef = NativeCRef<SDParamsData>;
+
+  public:
+    inline CELER_FUNCTION SensitiveDetectorView(SDParamsRef const& params);
+
+    inline auto CELER_FUNCTION detector_id(ImplVolumeId iv_id);
+
+  private:
+    SDParamsRef const& params_;
+};
+
+CELER_FUNCTION
+SensitiveDetectorView::SensitiveDetectorView(SDParamsRef const& params)
+    : params_(params)
+{
+    CELER_EXPECT(params_);
+}
+
+CELER_FUNCTION auto SensitiveDetectorView::detector_id(ImplVolumeId iv_id)
+{
+    return params_.detector[iv_id];
+}
+
+}  // namespace optical
+}  // namespace celeritas

src/celeritas/optical/surface/detail/PostBoundaryExecutor.hh

@@ -38,13 +38,28 @@ struct PostBoundaryExecutor
  */
 CELER_FUNCTION void PostBoundaryExecutor::operator()(CoreTrackView& track) const
 {
+    auto geo = track.geometry();
+    auto det_params = track.sensitive_detectors();
+    ImplVolumeId iv_id = geo.impl_volume_id();
+
+    DetectorId det_id = det_params.detector_id(iv_id);
+    if (det_id)
+    {
+        auto energy = track.particle().energy();
+        // TODO print energy when killing at SD
+        std::cout << "Detector hit in volume " << iv_id.get()
+                  << " on detector " << det_id.get() << " with energy "
+                  << energy.value() << std::endl;
+        track.sim().status(TrackStatus::killed);
+    }
+
     auto traverse = track.surface_physics().traversal();
     CELER_EXPECT(traverse.is_exiting());
 
     if (traverse.in_pre_volume())
     {
         // Re-entrant into the pre-volume
-        auto geo = track.geometry();
+        // auto geo = track.geometry();
         geo.cross_boundary();
         if (CELER_UNLIKELY(geo.failed()))
         {

src/celeritas/setup/Model.cc

@@ -14,6 +14,7 @@
 #include "geocel/VolumeParams.hh"
 #include "geocel/inp/Model.hh"
 #include "celeritas/geo/CoreGeoParams.hh"  // IWYU pragma: keep
+#include "celeritas/user/SDParams.hh"
 
 namespace celeritas
 {
@@ -118,6 +119,8 @@ ModelLoaded model(inp::Model const& m)
     result.surface
         = std::make_shared<SurfaceParams>(m.surfaces, *result.volume);
 
+    result.detector = std::make_shared<SDParams>(*result.geometry, m.detectors);
+
     return result;
 }
 

src/celeritas/setup/Model.hh

@@ -20,6 +20,7 @@ struct Model;
 
 class SurfaceParams;
 class VolumeParams;
+class SDParams;
 
 namespace setup
 {
@@ -30,6 +31,7 @@ struct ModelLoaded
     std::shared_ptr<CoreGeoParams> geometry;
     std::shared_ptr<SurfaceParams> surface;
     std::shared_ptr<VolumeParams> volume;
+    std::shared_ptr<SDParams> detector;
 };
 
 //---------------------------------------------------------------------------//

src/celeritas/setup/Problem.cc

@@ -313,6 +313,7 @@ auto build_optical_params(inp::OpticalPhysics const& optical,
     params.action_reg = std::make_shared<ActionRegistry>();
     params.gen_reg = std::make_shared<GeneratorRegistry>();
     params.max_streams = core.max_streams();
+    params.detectors = core.detector();
     {
         // Construct optical physics models
         optical::PhysicsParams::Input pp_inp;
@@ -459,6 +460,7 @@ ProblemLoaded problem(inp::Problem const& p, ImportData const& imported)
             }
             params.surface = std::make_shared<SurfaceParams>();
         }
+        params.detector = std::move(loaded_model.detector);
     }
 
     // Load materials