Prune primaries that start outside the world volume (#1624)

* Rename accumultaed counts in local transporter

* Discard tracks outside the world

* Summarize lost primaries at end of run

src/accel/LocalTransporter.cc

@@ -23,6 +23,7 @@
 
 #include "corecel/Config.hh"
 
+#include "corecel/cont/ArrayIO.hh"
 #include "corecel/cont/Span.hh"
 #include "corecel/io/Logger.hh"
 #include "corecel/sys/Device.hh"
@@ -96,6 +97,7 @@ LocalTransporter::LocalTransporter(SetupOptions const& options,
                       "offloading is disabled");
     particles_ = params.Params()->particle();
     CELER_ASSERT(particles_);
+    bbox_ = params.bbox();
 
     auto thread_id = get_geant_thread_id();
     CELER_VALIDATE(thread_id >= 0,
@@ -171,7 +173,7 @@ void LocalTransporter::InitializeEvent(int id)
     CELER_EXPECT(id >= 0);
 
     event_id_ = id_cast<UniqueEventId>(id);
-    ++accum_num_events_;
+    ++run_accum_.events;
 
     if (!(G4Threading::IsMultithreadedApplication()
           && G4MTRunManager::SeedOncePerCommunication()))
@@ -191,6 +193,23 @@ void LocalTransporter::Push(G4Track const& g4track)
 {
     CELER_EXPECT(*this);
 
+    if (Real3 pos = convert_from_geant(g4track.GetPosition(), 1);
+        !is_inside(bbox_, pos))
+    {
+        // Primary may have been created by a particle generator outside the
+        // geometry
+        double energy = g4track.GetKineticEnergy() / CLHEP::MeV;
+        CELER_LOG(debug)
+            << "Discarding track outside world: " << energy << " MeV from "
+            << g4track.GetDefinition()->GetParticleName() << " at " << pos
+            << " along "
+            << convert_from_geant(g4track.GetMomentumDirection(), 1);
+
+        buffer_accum_.lost_energy += energy;
+        ++buffer_accum_.lost_primaries;
+        return;
+    }
+
     Primary track;
 
     PDGNumber const pdg{g4track.GetDefinition()->GetPDGEncoding()};
@@ -221,7 +240,7 @@ void LocalTransporter::Push(G4Track const& g4track)
     track.event_id = EventId{0};
 
     buffer_.push_back(track);
-    buffer_energy_ += track.energy.value();
+    buffer_accum_.energy += track.energy.value();
     if (buffer_.size() >= auto_flush_)
     {
         /*!
@@ -268,9 +287,19 @@ void LocalTransporter::Flush()
     {
         CELER_LOG_LOCAL(debug)
             << "Transporting " << buffer_.size() << " tracks ("
-            << buffer_energy_ << " MeV cumulative kinetic energy) from event "
+            << buffer_accum_.energy
+            << " MeV cumulative kinetic energy) from event "
             << event_id_.unchecked_get() << " with Celeritas";
     }
+    if (buffer_accum_.lost_primaries > 0)
+    {
+        CELER_LOG_LOCAL(info)
+            << "Lost " << buffer_accum_.lost_energy
+            << " MeV cumulative kinetic energy from "
+            << buffer_accum_.lost_primaries
+            << " primaries that started outside the geometry in event "
+            << event_id_.unchecked_get();
+    }
 
     if (dump_primaries_)
     {
@@ -289,11 +318,12 @@ void LocalTransporter::Flush()
 
     // Copy buffered tracks to device and transport the first step
     auto track_counts = (*step_)(make_span(buffer_));
-    accum_num_steps_ += track_counts.active;
-    accum_num_primaries_ += buffer_.size();
+    run_accum_.steps += track_counts.active;
+    run_accum_.primaries += buffer_.size();
+    run_accum_.lost_primaries += buffer_accum_.lost_primaries;
 
     buffer_.clear();
-    buffer_energy_ = 0;
+    buffer_accum_ = {};
 
     size_type step_iters = 1;
 
@@ -306,7 +336,7 @@ void LocalTransporter::Flush()
                                       *step_);
 
         track_counts = (*step_)();
-        accum_num_steps_ += track_counts.active;
+        run_accum_.steps += track_counts.active;
         ++step_iters;
 
         CELER_VALIDATE_OR_KILL_ACTIVE(
@@ -328,10 +358,16 @@ void LocalTransporter::Finalize()
                    << "offloaded tracks (" << buffer_.size()
                    << " in buffer) were not flushed");
 
-    CELER_LOG_LOCAL(info) << "Finalizing Celeritas after " << accum_num_steps_
-                          << " steps from " << accum_num_primaries_
-                          << " offloaded tracks over " << accum_num_events_
+    CELER_LOG_LOCAL(info) << "Finalizing Celeritas after " << run_accum_.steps
+                          << " steps from " << run_accum_.primaries
+                          << " offloaded tracks over " << run_accum_.events
                           << " events";
+    if (run_accum_.lost_primaries > 0)
+    {
+        CELER_LOG_LOCAL(warning)
+            << "Lost a total of " << run_accum_.lost_primaries
+            << " primaries that started outside the world";
+    }
 
     if constexpr (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_GEANT4)
     {

src/accel/LocalTransporter.hh

@@ -13,6 +13,7 @@
 
 #include "corecel/Types.hh"
 #include "corecel/io/Logger.hh"
+#include "geocel/BoundingBox.hh"
 #include "celeritas/Types.hh"
 #include "celeritas/global/CoreParams.hh"
 #include "celeritas/global/Stepper.hh"
@@ -92,9 +93,32 @@ class LocalTransporter
     explicit operator bool() const { return static_cast<bool>(step_); }
 
   private:
+    //// TYPES ////
+
     using SPOffloadWriter = std::shared_ptr<detail::OffloadWriter>;
+    using BBox = BoundingBox<double>;
+
+    struct BufferAccum
+    {
+        double energy{0};  // MeV
+        double lost_energy{0};  // MeV
+        std::size_t lost_primaries{0};
+    };
+
+    struct RunAccum
+    {
+        std::size_t events{0};
+        std::size_t primaries{0};
+        std::size_t steps{0};
+        std::size_t lost_primaries{0};
+    };
+
+    //// DATA ////
 
     std::shared_ptr<ParticleParams const> particles_;
+    BBox bbox_;
+
+    // Thread-local data
     std::shared_ptr<StepperInterface> step_;
     std::vector<Primary> buffer_;
     std::shared_ptr<detail::HitProcessor> hit_processor_;
@@ -105,11 +129,9 @@ class LocalTransporter
 
     size_type auto_flush_{};
     size_type max_step_iters_{};
-    double buffer_energy_{0};
 
-    std::size_t accum_num_events_{0};
-    std::size_t accum_num_primaries_{0};
-    std::size_t accum_num_steps_{0};
+    BufferAccum buffer_accum_;
+    RunAccum run_accum_;
 
     // Shared across threads to write flushed particles
     SPOffloadWriter dump_primaries_;

src/accel/SharedParams.cc

@@ -20,11 +20,13 @@
 #include <G4Positron.hh>
 #include <G4RunManager.hh>
 #include <G4Threading.hh>
+#include <G4VisExtent.hh>
 
 #include "corecel/Config.hh"
 #include "corecel/Version.hh"
 
 #include "corecel/Assert.hh"
+#include "corecel/cont/ArrayIO.hh"
 #include "corecel/io/BuildOutput.hh"
 #include "corecel/io/Logger.hh"
 #include "corecel/io/OutputInterfaceAdapter.hh"
@@ -320,6 +322,24 @@ SharedParams::SharedParams(SetupOptions const& options)
     // Translate supported particles
     particles_ = build_g4_particles(params_->particle(), params_->physics());
 
+    // Create bounding box from navigator geometry
+    bbox_ = [] {
+        G4VPhysicalVolume const* world = GeantImporter::get_world_volume();
+        CELER_ASSERT(world);
+        G4LogicalVolume const* world_lv = world->GetLogicalVolume();
+        CELER_ASSERT(world_lv);
+        G4VSolid const* solid = world_lv->GetSolid();
+        CELER_ASSERT(solid);
+        G4VisExtent const& extent = solid->GetExtent();
+
+        BBox result{{extent.GetXmin(), extent.GetYmin(), extent.GetZmin()},
+                    {extent.GetXmax(), extent.GetYmax(), extent.GetZmax()}};
+        CELER_VALIDATE(result,
+                       << "world bounding box {" << result.lower() << ", "
+                       << result.upper() << "} is invalid");
+        return result;
+    }();
+
     // Save output filename (possibly empty if disabling output)
     output_filename_ = options.output_file;
 

src/accel/SharedParams.hh

@@ -11,6 +11,7 @@
 #include <vector>
 
 #include "corecel/Assert.hh"
+#include "geocel/BoundingBox.hh"
 
 class G4ParticleDefinition;
 
@@ -131,6 +132,7 @@ class SharedParams
     using SPOutputRegistry = std::shared_ptr<OutputRegistry>;
     using SPState = std::shared_ptr<CoreStateInterface>;
     using SPConstGeantGeoParams = std::shared_ptr<GeantGeoParams const>;
+    using BBox = BoundingBox<double>;
 
     //! Initialization status and integration mode
     Mode mode() const { return mode_; }
@@ -152,6 +154,9 @@ class SharedParams
 
     // Geant geometry wrapper, lazily created
     SPConstGeantGeoParams const& geant_geo_params() const;
+
+    // Geometry bounding box (CLHEP units)
+    BBox const& bbox() const { return bbox_; }
     //!@}
 
   private:
@@ -170,6 +175,7 @@ class SharedParams
     // Lazily created
     SPOutputRegistry output_reg_;
     SPConstGeantGeoParams geant_geo_;
+    BBox bbox_;
 
     //// HELPER FUNCTIONS ////