Allow detection in non-optical materials through an optical surface (#2308)

* Remove rindex and vestigial reflectivity

* Allow detector executor to detect hits on non-optical materials

* Apply feedback and add comments

* Delete old CI names

* Add build_sd_names virtual method to GeantTestBase

Prompt: "Add a virtual method to GeantTestBase that returns SetDetectors, then
passed to GeantSetup. Update the GeantTestBase::ImportSetup to include that
list. Then add the x/y/z detectors to the GDML file (auxtype="SensDet")
assuming the exising VolumeId is the order of the volume in that file."

Assisted-by: GitHub Copilot (claude-sonnet-4-5)

* Add note to GGParams

* Delete now-unnnecessary block and silence disabled physics

* Fix squirrel

Author: sethrj

src/celeritas/optical/CoreTrackView.hh

@@ -163,9 +163,6 @@ CoreTrackView::operator=(TrackInitializer const& init)
     // Initialize the surface state
     this->surface_physics().reset();
 
-    // Clear detector state data
-    states_.detectors.detector_hits[track_slot_id_].detector = {};
-
     return *this;
 }
 

src/celeritas/optical/action/detail/DetectorExecutor.hh

@@ -6,6 +6,7 @@
 //---------------------------------------------------------------------------//
 #pragma once
 
+#include "celeritas/Types.hh"
 #include "celeritas/optical/CoreTrackView.hh"
 #include "celeritas/optical/DetectorData.hh"
 
@@ -23,7 +24,14 @@ namespace detail
  * or is not in a detector region, an invalid hit is set in the corresponding
  * buffer track slot.
  *
+ * Detection can occur after a surface crossing \em or (for coupled EM
+ * problems) if an optical photon is emitted inside a detector region.
+ *
  * When a track generates a valid hit, it is killed (absorbed by the detector).
+ *
+ * \note This is called by \em all track slots, which is necessary to clear the
+ * outgoing detector hits. We could break this into a "reset" kernel that
+ * applies to all slots, and a "fill" kernel that applies only to valid slots.
  */
 struct DetectorExecutor
 {
@@ -43,41 +51,49 @@ CELER_FUNCTION void
 DetectorExecutor::operator()(CoreTrackView const& track) const
 {
     auto& hit = detector_state_.detector_hits[track.track_slot_id()];
+    // Clear the hit if inactive, errored, or not detected
+    hit.detector = {};
+
     auto sim = track.sim();
 
-    if (sim.status() == TrackStatus::alive)
+    if (!is_track_valid(sim.status()))
     {
-        auto const detectors = track.detectors();
-
-        auto geometry = track.geometry();
+        // Inactive or errored
+        return;
+    }
+    if (track.surface_physics().is_crossing_boundary())
+    {
+        // Boundary crossing not yet completed, so not yet detected!
+        return;
+    }
 
-        auto const volume_id = geometry.volume_id();
-        auto const detector_id = detectors.detector_id(volume_id);
+    auto geometry = track.geometry();
+    if (geometry.is_outside())
+    {
+        // Killed by leaving geometry; no detection
+        return;
+    }
 
-        if (detector_id)
-        {
-            // Score a valid hit
-            hit = DetectorHit{detector_id,
-                              track.sim().primary_id(),
-                              track.particle().energy(),
-                              sim.time(),
-                              geometry.pos(),
-                              geometry.volume_instance_id()};
+    auto const detectors = track.detectors();
+    auto const volume_id = geometry.volume_id();
+    auto const detector_id = detectors.detector_id(volume_id);
 
-            // Kill the track
-            sim.status(TrackStatus::killed);
-        }
-        else
-        {
-            // Mark that the track is not in a detector
-            hit.detector = {};
-        }
-    }
-    else
+    if (!detector_id)
     {
-        // Ensure killed, inactive, and errored tracks don't contribute to hits
-        hit.detector = {};
+        // Not in a detector
+        return;
     }
+
+    // Score a valid hit
+    hit.detector = detector_id;
+    hit.primary = sim.primary_id();
+    hit.energy = track.particle().energy();
+    hit.time = sim.time();
+    hit.position = geometry.pos();
+    hit.volume_instance = geometry.volume_instance_id();
+
+    // Kill the track
+    sim.status(TrackStatus::killed);
 }
 
 //---------------------------------------------------------------------------//

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

@@ -81,15 +81,7 @@ CELER_FUNCTION void InitBoundaryExecutor::operator()(CoreTrackView& track) const
         = select_surface(track.surface(), geo.volume_instance_id());
     if (!oriented_surface)
     {
-        // Kill the track if the post-volume doesn't have a valid optical
-        // material and there's no surface
-        if (!post_volume_material)
-        {
-            track.sim().status(TrackStatus::killed);
-            return;
-        }
-
-        // Use default surface data
+        // Use default surface properties: typically dielectric-dielectric
         oriented_surface.surface = surface_physics.scalars().default_surface;
         oriented_surface.orientation = SubsurfaceDirection::forward;
     }
@@ -102,12 +94,16 @@ CELER_FUNCTION void InitBoundaryExecutor::operator()(CoreTrackView& track) const
         global_normal = -global_normal;
     }
 
-    surface_physics
-        = SurfacePhysicsTrackView::Initializer{oriented_surface.surface,
-                                               oriented_surface.orientation,
-                                               global_normal,
-                                               pre_volume_material,
-                                               post_volume_material};
+    surface_physics = [&] {
+        SurfacePhysicsTrackView::Initializer init;
+        init.surface = oriented_surface.surface;
+        init.orientation = oriented_surface.orientation;
+        init.global_normal = global_normal;
+        init.pre_volume_material = pre_volume_material;
+        // Note that post-volume material may be null
+        init.post_volume_material = post_volume_material;
+        return init;
+    }();
 
     CELER_ASSERT(
         is_entering_surface(geo.dir(), surface_physics.global_normal()));

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

@@ -62,7 +62,8 @@ CELER_FUNCTION void PostBoundaryExecutor::operator()(CoreTrackView& track) const
 
     if (!track.material_record().material_id())
     {
-        // Kill track if it enters an invalid optical material
+        // Kill track if it enters an invalid optical material after crossing
+        // through a custom physics surface
         track.sim().status(TrackStatus::killed);
     }
 

src/geocel/GeantGeoParams.cc

@@ -720,7 +720,8 @@ GeantGeoParams::from_gdml(std::string const& filename)
     // Create geo params
     auto result
         = std::make_shared<GeantGeoParams>(loaded.world, Ownership::value);
-    // Set detectors (hack)
+    // We own constructed detectors (note that these live only on the main
+    // thread and are not suitable for G4 MT: use DetectorConstruction instead)
     result->built_detectors_ = std::move(built_detectors);
 
     // Save for use outside in Celeritas

test/celeritas/GeantTestBase.cc

@@ -42,6 +42,7 @@ struct GeantTestBase::ImportSetup
     std::shared_ptr<GeantGeoParams> geo;
     GeantPhysicsOptions options{};
     GeantImportDataSelection selection{};
+    GeantSetup::SetString sd_names{};
     ImportData imported;
     ScopedGeantExceptionHandler scoped_exceptions;
 };
@@ -64,21 +65,6 @@ bool GeantTestBase::is_ci_build()
     return clhep >= Version{2, 4, 6} && clhep < Version{2, 5}
            && g4 >= Version{11, 3} && g4 < Version{11, 4};
 }
-
-//---------------------------------------------------------------------------//
-//! Whether Geant4 dependencies match those on Wildstyle
-bool GeantTestBase::is_wildstyle_build()
-{
-    return GeantTestBase::is_ci_build();
-}
-
-//---------------------------------------------------------------------------//
-//! Whether Geant4 dependencies match those on Summit
-bool GeantTestBase::is_summit_build()
-{
-    return GeantTestBase::is_ci_build();
-}
-
 //---------------------------------------------------------------------------//
 // PROTECTED MEMBER FUNCTIONS
 //---------------------------------------------------------------------------//
@@ -146,6 +132,7 @@ auto GeantTestBase::load(std::string const& filename) const
 {
     GeantPhysicsOptions opts = this->build_geant_options();
     GeantImportDataSelection sel = this->build_import_data_selection();
+    GeantSetup::SetString sd_names = this->build_sd_names();
 
     using PersistentImportSetup = PersistentSP<GeantTestBase::ImportSetup>;
 
@@ -158,7 +145,9 @@ auto GeantTestBase::load(std::string const& filename) const
                        << "load was called before build_geant_geo");
         i = std::make_shared<ImportSetup>();
         i->options = opts;
-        i->import = std::make_unique<GeantImporter>(GeantSetup{filename, opts});
+        i->sd_names = sd_names;
+        i->import = std::make_unique<GeantImporter>(
+            GeantSetup{filename, opts, std::move(sd_names)});
         i->geo = i->import->geo_params();
         CELER_ASSERT(i->geo);
         CELER_ASSERT(!celeritas::global_geant_geo().expired());
@@ -181,6 +170,9 @@ auto GeantTestBase::load(std::string const& filename) const
         CELER_VALIDATE(opts == i->options,
                        << "cannot change physics options after setup "
                        << explanation);
+        CELER_VALIDATE(sd_names == i->sd_names,
+                       << "cannot change sensitive detector names after setup "
+                       << explanation);
 
         if (sel == i->selection)
         {
@@ -206,6 +198,12 @@ auto GeantTestBase::load(std::string const& filename) const
     return *i;
 }
 
+//---------------------------------------------------------------------------//
+auto GeantTestBase::build_sd_names() const -> GeantSetup::SetString
+{
+    return {};
+}
+
 //---------------------------------------------------------------------------//
 GeantImportDataSelection GeantTestBase::build_import_data_selection() const
 {

test/celeritas/GeantTestBase.hh

@@ -9,6 +9,7 @@
 #include <iosfwd>
 
 #include "celeritas/Types.hh"
+#include "celeritas/ext/GeantSetup.hh"
 
 #include "ImportedDataTestBase.hh"
 
@@ -17,7 +18,6 @@ class G4VPhysicalVolume;
 namespace celeritas
 {
 //---------------------------------------------------------------------------//
-struct GeantPhysicsOptions;
 struct GeantImportDataSelection;
 
 namespace test
@@ -34,12 +34,11 @@ class GeantTestBase : public ImportedDataTestBase
     //!@{
     //! Whether the Geant4 configuration match a certain machine
     static bool is_ci_build();
-    static bool is_wildstyle_build();
-    static bool is_summit_build();
     //!@}
 
   protected:
     virtual GeantPhysicsOptions build_geant_options() const;
+    virtual GeantSetup::SetString build_sd_names() const;
 
     SPConstTrackInit build_init() override;
     SPConstAction build_along_step() override;

test/celeritas/ext/GeantImporter.test.cc

@@ -457,8 +457,8 @@ TEST_F(DuneCryostat, optical_gen)
     ASSERT_TRUE(scint);
     EXPECT_EQ(1, scint->materials.size());
 
-    ASSERT_TRUE(scint->materials.count(OptMatId{1}));
-    auto const& m = scint->materials.at(OptMatId{1});
+    ASSERT_TRUE(scint->materials.count(OptMatId{0}));
+    auto const& m = scint->materials.at(OptMatId{0});
     ASSERT_EQ(2, m.components.size());
     EXPECT_SOFT_EQ(50000 * 0.8, m.components[0].yield);
     EXPECT_SOFT_EQ(50000 * 0.2, m.components[1].yield);

test/celeritas/optical/Detector.test.cc

@@ -56,6 +56,10 @@ class DetectorTest : public ::celeritas::test::GeantTestBase
     {
         auto result = GeantTestBase::build_geant_options();
         result.optical.emplace();
+        result.optical->scintillation = std::nullopt;
+        result.optical->mie_scattering = false;
+        result.optical->wavelength_shifting = std::nullopt;
+        result.optical->wavelength_shifting2 = std::nullopt;
         return result;
     }
 
@@ -66,6 +70,11 @@ class DetectorTest : public ::celeritas::test::GeantTestBase
         return result;
     }
 
+    GeantSetup::SetString build_sd_names() const override
+    {
+        return {"x-detectors", "y-detectors", "z-detectors"};
+    }
+
     std::vector<IMC> select_optical_models() const override
     {
         return {IMC::absorption};
@@ -87,23 +96,6 @@ class DetectorTest : public ::celeritas::test::GeantTestBase
             this->optical_action_reg().get(), input);
     }
 
-    SPConstDetectors detector() override
-    {
-        if (!detector_)
-        {
-            inp::Detectors input{{
-                {"y-detectors", {VolumeId{1}, VolumeId{2}}},
-                {"x-detectors", {VolumeId{3}, VolumeId{4}}},
-                {"z-detectors", {VolumeId{5}, VolumeId{6}}},
-            }};
-
-            detector_ = std::make_shared<DetectorParams>(std::move(input),
-                                                         *this->volume());
-        }
-
-        return detector_;
-    }
-
     inp::OpticalDetector build_optical_detector_input() override
     {
         return detector_input_;

test/celeritas/optical/Generator.test.cc

@@ -267,7 +267,7 @@ TEST_F(DuneGeneratorTest, offload)
     gdd.primary = PrimaryId{123};
     gdd.step_length = 2.0 * units::centimeter;
     gdd.charge = units::ElementaryCharge{-1};
-    gdd.material = OptMatId{1};  // Should be LAr
+    gdd.material = OptMatId{0};  // Should be LAr
     gdd.continuous_edep_fraction = 1.0;
     gdd.points[StepPoint::pre] = {units::LightSpeed(0.7),
                                   1e-9 * units::nanosecond,