Automatically calculate number of events in apps (#1072)

* Automatically calculate number of events in apps

* Address review feedback

* Fix counting of events from HepMC3 file

- `Reader::skip()` was added in HepMC3 version 3.2.0: use `read_event()` to
  count with older versions
- Fix count and validate for zero events

app/celer-g4/ActionInitialization.cc

@@ -36,21 +36,7 @@ ActionInitialization::ActionInitialization(SPParams params)
     // Create Geant4 diagnostics to be shared across worker threads
     diagnostics_ = std::make_shared<GeantDiagnostics>();
 
-    auto const& inp = GlobalSetup::Instance()->input();
-    CELER_VALIDATE(inp.primary_options || !inp.event_file.empty(),
-                   << "no event input file nor primary options were "
-                      "specified");
-
-    if (!inp.event_file.empty())
-    {
-        hepmc_gen_ = std::make_shared<HepMC3PrimaryGenerator>(inp.event_file);
-        num_events_ = hepmc_gen_->NumEvents();
-    }
-    else
-    {
-        num_events_ = inp.primary_options.num_events;
-    }
-
+    num_events_ = GlobalSetup::Instance()->setup_options().max_num_events;
     CELER_ENSURE(num_events_ > 0);
 }
 
@@ -87,9 +73,9 @@ void ActionInitialization::Build() const
     CELER_LOG_LOCAL(status) << "Constructing user action";
 
     // Primary generator emits source particles
-    if (hepmc_gen_)
+    if (auto hepmc_gen = GlobalSetup::Instance()->hepmc_gen())
     {
-        this->SetUserAction(new HepMC3PrimaryGeneratorAction(hepmc_gen_));
+        this->SetUserAction(new HepMC3PrimaryGeneratorAction(hepmc_gen));
     }
     else
     {

app/celer-g4/ActionInitialization.hh

@@ -16,7 +16,6 @@
 
 namespace celeritas
 {
-class HepMC3PrimaryGenerator;
 namespace app
 {
 //---------------------------------------------------------------------------//
@@ -43,7 +42,6 @@ class ActionInitialization final : public G4VUserActionInitialization
   private:
     SPParams params_;
     SPDiagnostics diagnostics_;
-    std::shared_ptr<HepMC3PrimaryGenerator> hepmc_gen_;
     int num_events_{0};
     mutable bool init_shared_;
 };

app/celer-g4/GlobalSetup.cc

@@ -19,6 +19,7 @@
 #include "celeritas/ext/RootFileManager.hh"
 #include "celeritas/field/RZMapFieldInput.hh"
 #include "accel/ExceptionConverter.hh"
+#include "accel/HepMC3PrimaryGenerator.hh"
 #include "accel/SetupOptionsMessenger.hh"
 
 #include "HepMC3PrimaryGeneratorAction.hh"
@@ -151,7 +152,21 @@ void GlobalSetup::ReadInput(std::string const& filename)
 
         // Apply Celeritas \c SetupOptions commands
         options_->max_num_tracks = input_.num_track_slots;
-        options_->max_num_events = input_.max_events;
+        options_->max_num_events = [this] {
+            CELER_VALIDATE(input_.primary_options || !input_.event_file.empty(),
+                           << "no event input file nor primary options were "
+                              "specified");
+            if (!input_.event_file.empty())
+            {
+                hepmc_gen_ = std::make_shared<HepMC3PrimaryGenerator>(
+                    input_.event_file);
+                return static_cast<size_type>(hepmc_gen_->NumEvents());
+            }
+            else
+            {
+                return input_.primary_options.num_events;
+            }
+        }();
         options_->max_steps = input_.max_steps;
         options_->initializer_capacity = input_.initializer_capacity;
         options_->secondary_stack_factor = input_.secondary_stack_factor;

app/celer-g4/GlobalSetup.hh

@@ -23,6 +23,7 @@ class G4GenericMessenger;
 
 namespace celeritas
 {
+class HepMC3PrimaryGenerator;
 namespace app
 {
 //---------------------------------------------------------------------------//
@@ -31,6 +32,12 @@ namespace app
  */
 class GlobalSetup
 {
+  public:
+    //!@{
+    //! \name Type aliases
+    using SPPrimaryGenerator = std::shared_ptr<HepMC3PrimaryGenerator>;
+    //!@}
+
   public:
     // Return non-owning pointer to a singleton
     static GlobalSetup* Instance();
@@ -91,13 +98,17 @@ class GlobalSetup
     //! Whether ROOT I/O for SDs is enabled
     bool root_sd_io() const { return root_sd_io_; }
 
+    //! Get HepMC3 primary generator
+    SPPrimaryGenerator hepmc_gen() const { return hepmc_gen_; }
+
   private:
     // Private constructor since we're a singleton
     GlobalSetup();
     ~GlobalSetup();
 
     // Data
     std::shared_ptr<SetupOptions> options_;
+    std::shared_ptr<HepMC3PrimaryGenerator> hepmc_gen_;
     RunInput input_;
     Stopwatch get_setup_time_;
     bool root_sd_io_{false};

app/celer-g4/RunInput.cc

@@ -50,7 +50,7 @@ RunInput::operator bool() const
     return !geometry_file.empty() && (primary_options || !event_file.empty())
            && physics_list < PhysicsListSelection::size_
            && (field == no_field() || field_options)
-           && ((num_track_slots > 0 && max_events > 0 && max_steps > 0
+           && ((num_track_slots > 0 && max_steps > 0
                 && initializer_capacity > 0 && secondary_stack_factor > 0)
                || SharedParams::CeleritasDisabled())
            && (step_diagnostic_bins > 0 || !step_diagnostic);

app/celer-g4/RunInput.hh

@@ -64,7 +64,6 @@ struct RunInput
 
     // Control
     size_type num_track_slots{};
-    size_type max_events{};
     size_type max_steps{unspecified};
     size_type initializer_capacity{};
     real_type secondary_stack_factor{};

app/celer-g4/RunInputIO.json.cc

@@ -67,7 +67,6 @@ void from_json(nlohmann::json const& j, RunInput& v)
     RI_LOAD_OPTION(primary_options);
 
     RI_LOAD_OPTION(num_track_slots);
-    RI_LOAD_OPTION(max_events);
     RI_LOAD_OPTION(max_steps);
     RI_LOAD_OPTION(initializer_capacity);
     RI_LOAD_OPTION(secondary_stack_factor);
@@ -156,7 +155,6 @@ void to_json(nlohmann::json& j, RunInput const& v)
     }
 
     RI_SAVE(num_track_slots);
-    RI_SAVE(max_events);
     RI_SAVE_OPTION(max_steps);
     RI_SAVE(initializer_capacity);
     RI_SAVE(secondary_stack_factor);

app/celer-g4/test-harness.py

@@ -71,7 +71,6 @@ def strtobool(text):
     "output_file": out_file,
     "offload_output_file": offload_file,
     "num_track_slots": max_tracks,
-    "max_events": 1024,
     "initializer_capacity": init_capacity,
     "secondary_stack_factor": 2,
     "physics_list": "ftfp_bert",

app/celer-sim/Runner.cc

@@ -84,7 +84,7 @@ namespace
  * value, the number of threads for each nested parallel region will be set to
  * that value.
  */
-size_type calc_num_streams(RunnerInput const& inp)
+size_type calc_num_streams(RunnerInput const& inp, size_type num_events)
 {
     size_type num_threads = 1;
 #if CELERITAS_USE_OPENMP
@@ -102,7 +102,7 @@ size_type calc_num_streams(RunnerInput const& inp)
     CELER_DISCARD(inp);
 #endif
     // Don't create more streams than events
-    return std::min(num_threads, inp.max_events);
+    return std::min(num_threads, num_events);
 }
 
 //---------------------------------------------------------------------------//
@@ -123,7 +123,6 @@ Runner::Runner(RunnerInput const& inp, SPOutputRegistry output)
     this->build_diagnostics(inp);
     this->build_step_collectors(inp);
     this->build_transporter_input(inp);
-    this->build_events(inp);
     use_device_ = inp.use_device;
 
     if (root_manager_)
@@ -368,28 +367,24 @@ void Runner::build_core_params(RunnerInput const& inp,
     // Construct simulation params
     params.sim = SimParams::from_import(imported, params.particle);
 
+    // Get the total number of events
+    auto num_events = this->build_events(inp, params.particle);
+
     // Store the number of simultaneous threads/tasks per process
-    params.max_streams = calc_num_streams(inp);
+    params.max_streams = calc_num_streams(inp, num_events);
     CELER_VALIDATE(inp.mctruth_file.empty() || params.max_streams == 1,
                    << "cannot output MC truth with multiple "
                       "streams ("
                    << params.max_streams << " requested)");
 
     // Construct track initialization params
-    params.init = [&inp, &params] {
+    params.init = [&inp, &params, num_events] {
         CELER_VALIDATE(inp.initializer_capacity > 0,
                        << "nonpositive initializer_capacity="
                        << inp.initializer_capacity);
-        CELER_VALIDATE(inp.max_events > 0,
-                       << "nonpositive max_events=" << inp.max_events);
-        CELER_VALIDATE(!inp.primary_options
-                           || inp.max_events >= inp.primary_options.num_events,
-                       << "max_events=" << inp.max_events
-                       << " cannot be less than num_events="
-                       << inp.primary_options.num_events);
         TrackInitParams::Input input;
         input.capacity = ceil_div(inp.initializer_capacity, params.max_streams);
-        input.max_events = inp.max_events;
+        input.max_events = num_events;
         input.track_order = inp.track_order;
         return std::make_shared<TrackInitParams>(std::move(input));
     }();
@@ -421,8 +416,11 @@ void Runner::build_transporter_input(RunnerInput const& inp)
 //---------------------------------------------------------------------------//
 /*!
  * Read events from a file or build using a primary generator.
+ *
+ * This returns the total number of events.
  */
-void Runner::build_events(RunnerInput const& inp)
+size_type
+Runner::build_events(RunnerInput const& inp, SPConstParticles particles)
 {
     ScopedMem record_mem("Runner.build_events");
 
@@ -447,21 +445,22 @@ void Runner::build_events(RunnerInput const& inp)
             }
             event = generate();
         }
+        return generate.num_events();
     };
 
     if (inp.primary_options)
     {
-        read_events(PrimaryGenerator::from_options(core_params_->particle(),
-                                                   inp.primary_options));
+        return read_events(
+            PrimaryGenerator::from_options(particles, inp.primary_options));
     }
     else if (ends_with(inp.event_file, ".root"))
     {
-        read_events(RootEventReader(inp.event_file, core_params_->particle()));
+        return read_events(RootEventReader(inp.event_file, particles));
     }
     else
     {
         // Assume filename is one of the HepMC3-supported extensions
-        read_events(EventReader(inp.event_file, core_params_->particle()));
+        return read_events(EventReader(inp.event_file, particles));
     }
 }
 

app/celer-sim/Runner.hh

@@ -22,6 +22,7 @@ namespace celeritas
 {
 class CoreParams;
 class OutputRegistry;
+class ParticleParams;
 class RootFileManager;
 class StepCollector;
 }  // namespace celeritas
@@ -77,6 +78,7 @@ class Runner
     //// TYPES ////
 
     using UPTransporterBase = std::unique_ptr<TransporterBase>;
+    using SPConstParticles = std::shared_ptr<ParticleParams const>;
     using VecPrimary = std::vector<Primary>;
     using VecEvent = std::vector<VecPrimary>;
 
@@ -99,8 +101,7 @@ class Runner
     void build_step_collectors(RunnerInput const&);
     void build_diagnostics(RunnerInput const&);
     void build_transporter_input(RunnerInput const&);
-    void build_events(RunnerInput const&);
-    int get_num_streams(RunnerInput const&);
+    size_type build_events(RunnerInput const&, SPConstParticles);
     TransporterBase& get_transporter(StreamId);
     TransporterBase const* get_transporter_ptr(StreamId) const;
 };

app/celer-sim/RunnerInput.hh

@@ -68,7 +68,6 @@ struct RunnerInput
     size_type num_track_slots{};  //!< Divided among streams
     size_type max_steps{unspecified};
     size_type initializer_capacity{};  //!< Divided among streams
-    size_type max_events{};
     real_type secondary_stack_factor{};
     bool use_device{};
     bool sync{};
@@ -99,8 +98,7 @@ struct RunnerInput
         return !geometry_file.empty()
                && (primary_options || !event_file.empty())
                && num_track_slots > 0 && max_steps > 0
-               && initializer_capacity > 0 && max_events > 0
-               && secondary_stack_factor > 0
+               && initializer_capacity > 0 && secondary_stack_factor > 0
                && (step_diagnostic_bins > 0 || !step_diagnostic)
                && (field == no_field() || field_options);
     }

app/celer-sim/RunnerInputIO.json.cc

@@ -98,7 +98,6 @@ void from_json(nlohmann::json const& j, RunnerInput& v)
     LDIO_LOAD_OPTION(num_track_slots);
     LDIO_LOAD_OPTION(max_steps);
     LDIO_LOAD_REQUIRED(initializer_capacity);
-    LDIO_LOAD_REQUIRED(max_events);
     LDIO_LOAD_REQUIRED(secondary_stack_factor);
     LDIO_LOAD_REQUIRED(use_device);
     LDIO_LOAD_OPTION(sync);
@@ -173,7 +172,6 @@ void to_json(nlohmann::json& j, RunnerInput const& v)
     LDIO_SAVE(num_track_slots);
     LDIO_SAVE_OPTION(max_steps);
     LDIO_SAVE(initializer_capacity);
-    LDIO_SAVE(max_events);
     LDIO_SAVE(secondary_stack_factor);
     LDIO_SAVE(use_device);
     LDIO_SAVE(sync);

app/celer-sim/simple-driver.py

@@ -91,7 +91,6 @@ def strtobool(text):
     'num_track_slots': num_tracks,
     'max_steps': max_steps,
     'initializer_capacity': 100 * max([num_tracks, num_primaries]),
-    'max_events': 1000,
     'secondary_stack_factor': 3,
     'action_diagnostic': True,
     'step_diagnostic': True,

src/accel/HepMC3PrimaryGenerator.cc

@@ -115,11 +115,24 @@ HepMC3PrimaryGenerator::HepMC3PrimaryGenerator(std::string const& filename)
         SPReader temp_reader = open_hepmc3(filename);
         CELER_ASSERT(temp_reader);
         size_type result = 0;
-        while (!temp_reader->failed())
+#if HEPMC3_VERSION_CODE < 3002000
+        HepMC3::GenEvent evt;
+        temp_reader->read_event(evt);
+#else
+        temp_reader->skip(0);
+#endif
+        CELER_VALIDATE(!temp_reader->failed(),
+                       << "event file '" << filename
+                       << "' did not contain any events");
+        do
         {
-            temp_reader->skip(1);
             result++;
-        }
+#if HEPMC3_VERSION_CODE < 3002000
+            temp_reader->read_event(evt);
+#else
+            temp_reader->skip(1);
+#endif
+        } while (!temp_reader->failed());
         CELER_LOG(debug) << "HepMC3 file has " << result << " events";
         return result;
     }();

src/celeritas/io/EventIOInterface.hh

@@ -10,6 +10,7 @@
 #include <vector>
 
 #include "corecel/Macros.hh"
+#include "corecel/Types.hh"
 
 namespace celeritas
 {
@@ -60,6 +61,9 @@ class EventReaderInterface
     //! Read all primaries from a single event
     virtual result_type operator()() = 0;
 
+    //! Get total number of events
+    virtual size_type num_events() const = 0;
+
   protected:
     //!@{
     //! Allow construction and assignment only through daughter classes