Change PrimaryGeneratorAction to be a wrapper (#1593)

* Add error checking to particle params
* Make accessible helper functions to convert particle/volume to ImportX
* Add helper to build particle params from import
* Add seed function to primary generator
* Make PG action into a "primary converter" wrapper function

app/celer-g4/ActionInitialization.cc

@@ -7,6 +7,7 @@
 #include "ActionInitialization.hh"
 
 #include "corecel/io/Logger.hh"
+#include "celeritas/phys/PrimaryGeneratorOptions.hh"
 #include "accel/ExceptionConverter.hh"
 #include "accel/HepMC3PrimaryGenerator.hh"
 #include "accel/LocalTransporter.hh"

app/celer-g4/PGPrimaryGeneratorAction.cc

@@ -13,34 +13,58 @@
 #include "corecel/Macros.hh"
 #include "geocel/GeantUtils.hh"
 #include "geocel/g4/Convert.hh"
-#include "celeritas/phys/PrimaryGeneratorOptions.hh"
+#include "celeritas/ext/GeantImporter.hh"
+#include "celeritas/ext/GeantUnits.hh"
+#include "celeritas/phys/ParticleParams.hh"
+#include "celeritas/phys/Primary.hh"
 
 namespace celeritas
 {
 namespace app
 {
+namespace
+{
+//---------------------------------------------------------------------------//
+auto make_particles(std::vector<PDGNumber> const& all_pdg)
+{
+    CELER_VALIDATE(!all_pdg.empty(),
+                   << "primary generator has no input particles");
+
+    auto* par_table = G4ParticleTable::GetParticleTable();
+    CELER_ASSERT(par_table);
+
+    // Find and convert Geant4 particles
+    ParticleParams::Input inp;
+    for (auto pdg : all_pdg)
+    {
+        CELER_EXPECT(pdg);
+        auto* p = par_table->FindParticle(pdg.get());
+        CELER_VALIDATE(
+            p, << "particle with PDG " << pdg.get() << " is not loaded");
+        inp.push_back(
+            ParticleParams::ParticleInput::from_import(import_particle(*p)));
+    }
+
+    return std::make_shared<ParticleParams>(std::move(inp));
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace
+
 //---------------------------------------------------------------------------//
 /*!
  * Construct primary action.
  */
-PGPrimaryGeneratorAction::PGPrimaryGeneratorAction(
-    PrimaryGeneratorOptions const& options)
+PGPrimaryGeneratorAction::PGPrimaryGeneratorAction(Input const& i)
+    : particle_params_{make_particles(i.pdg)}
+    , generate_primaries_{PrimaryGenerator::from_options(particle_params_, i)}
 {
-    CELER_EXPECT(options);
-
     // Generate one particle at each call to \c GeneratePrimaryVertex()
     gun_.SetNumberOfParticles(1);
 
-    seed_ = options.seed;
-    num_events_ = options.num_events;
-    primaries_per_event_ = options.primaries_per_event;
-    sample_energy_ = make_energy_sampler(options.energy);
-    sample_pos_ = make_position_sampler(options.position);
-    sample_dir_ = make_direction_sampler(options.direction);
-
-    // Set the particle definitions
-    particle_def_.reserve(options.pdg.size());
-    for (auto const& pdg : options.pdg)
+    // Save the particle definitions corresponding to particle IDs
+    particle_def_.reserve(i.pdg.size());
+    for (auto const& pdg : i.pdg)
     {
         particle_def_.push_back(
             G4ParticleTable::GetParticleTable()->FindParticle(pdg.get()));
@@ -56,28 +80,23 @@ void PGPrimaryGeneratorAction::GeneratePrimaries(G4Event* event)
     CELER_EXPECT(event);
     CELER_EXPECT(event->GetEventID() >= 0);
 
-    size_type event_id = event->GetEventID();
-    if (event_id >= num_events_)
-    {
-        return;
-    }
-
     // Seed with an independent value for each event. Since Geant4 schedules
     // events dynamically, the same event ID may not be mapped to the same
     // thread across multiple runs. For reproducibility, Geant4 reseeds each
     // worker thread's RNG at the start of each event with a seed calculated
     // from the event ID.
-    rng_.seed(seed_ + event_id);
+    generate_primaries_.seed(id_cast<UniqueEventId>(event->GetEventID()));
+
+    auto primaries = generate_primaries_();
 
-    for (size_type i = 0; i < primaries_per_event_; ++i)
+    for (auto const& p : primaries)
     {
-        gun_.SetParticleDefinition(particle_def_[i % particle_def_.size()]);
-        gun_.SetParticlePosition(
-            convert_to_geant(sample_pos_(rng_), clhep_length));
-        gun_.SetParticleMomentumDirection(
-            convert_to_geant(sample_dir_(rng_), 1));
-        gun_.SetParticleEnergy(
-            convert_to_geant(sample_energy_(rng_), CLHEP::MeV));
+        CELER_ASSERT(p.particle_id < particle_def_.size());
+        gun_.SetParticleDefinition(
+            particle_def_[p.particle_id.unchecked_get()]);
+        gun_.SetParticlePosition(convert_to_geant(p.position, clhep_length));
+        gun_.SetParticleMomentumDirection(convert_to_geant(p.direction, 1));
+        gun_.SetParticleEnergy(convert_to_geant(p.energy, CLHEP::MeV));
         gun_.GeneratePrimaryVertex(event);
 
         if (CELERITAS_DEBUG)
@@ -88,7 +107,7 @@ void PGPrimaryGeneratorAction::GeneratePrimaries(G4Event* event)
     }
 
     CELER_ENSURE(event->GetNumberOfPrimaryVertex()
-                 == static_cast<int>(primaries_per_event_));
+                 == static_cast<int>(primaries.size()));
 }
 
 //---------------------------------------------------------------------------//

app/celer-g4/PGPrimaryGeneratorAction.hh

@@ -6,16 +6,14 @@
 //---------------------------------------------------------------------------//
 #pragma once
 
-#include <functional>
 #include <G4Event.hh>
 #include <G4ParticleDefinition.hh>
 #include <G4ParticleGun.hh>
 #include <G4VUserPrimaryGeneratorAction.hh>
 
 #include "corecel/Types.hh"
 #include "corecel/cont/Array.hh"
-#include "geocel/Types.hh"
-#include "celeritas/phys/PrimaryGeneratorOptions.hh"
+#include "celeritas/phys/PrimaryGenerator.hh"
 
 namespace celeritas
 {
@@ -36,28 +34,23 @@ class PGPrimaryGeneratorAction final : public G4VUserPrimaryGeneratorAction
   public:
     //!@{
     //! \name Type aliases
-    using EnergySampler = std::function<real_type(PrimaryGeneratorEngine&)>;
-    using PositionSampler = std::function<Real3(PrimaryGeneratorEngine&)>;
-    using DirectionSampler = std::function<Real3(PrimaryGeneratorEngine&)>;
+    using Input = PrimaryGeneratorOptions;
     //!@}
 
   public:
     // Construct primary action
-    explicit PGPrimaryGeneratorAction(PrimaryGeneratorOptions const& opts);
+    explicit PGPrimaryGeneratorAction(Input const& opts);
 
     // Generate events
     void GeneratePrimaries(G4Event* event) final;
 
   private:
+    using GeneratorImpl = celeritas::PrimaryGenerator;
+
+    GeneratorImpl::SPConstParticles particle_params_;
+    GeneratorImpl generate_primaries_;
     G4ParticleGun gun_;
-    PrimaryGeneratorEngine rng_;
-    size_type num_events_{};
-    size_type primaries_per_event_{};
     std::vector<G4ParticleDefinition*> particle_def_;
-    EnergySampler sample_energy_;
-    PositionSampler sample_pos_;
-    DirectionSampler sample_dir_;
-    size_type seed_{0};
 };
 
 //---------------------------------------------------------------------------//

src/celeritas/ext/GeantImporter.cc

@@ -329,36 +329,18 @@ import_particles(GeantImporter::DataSelection::Flags particle_flags)
 
     std::vector<ImportParticle> particles;
 
-    double const time_scale = native_value_from_clhep(ImportUnits::time);
-
     ParticleFilter include_particle{particle_flags};
     while (particle_iterator())
     {
-        G4ParticleDefinition const& g4_particle_def
-            = *(particle_iterator.value());
+        G4ParticleDefinition const& p = *(particle_iterator.value());
 
-        PDGNumber pdg{g4_particle_def.GetPDGEncoding()};
+        PDGNumber pdg{p.GetPDGEncoding()};
         if (!include_particle(pdg))
         {
             continue;
         }
 
-        ImportParticle particle;
-        particle.name = g4_particle_def.GetParticleName();
-        particle.pdg = pdg.unchecked_get();
-        particle.mass = g4_particle_def.GetPDGMass();
-        particle.charge = g4_particle_def.GetPDGCharge();
-        particle.spin = g4_particle_def.GetPDGSpin();
-        particle.lifetime = g4_particle_def.GetPDGLifeTime();
-        particle.is_stable = g4_particle_def.GetPDGStable();
-
-        if (!particle.is_stable)
-        {
-            // Convert lifetime of unstable particles to seconds
-            particle.lifetime *= time_scale;
-        }
-
-        particles.push_back(particle);
+        particles.push_back(import_particle(p));
     }
     CELER_LOG(debug) << "Loaded " << particles.size() << " particles";
     CELER_ENSURE(!particles.empty());
@@ -1208,7 +1190,7 @@ ImportData GeantImporter::operator()(DataSelection const& selected)
             }
         }
         imported.regions = import_regions();
-        imported.volumes = this->import_volumes(selected.unique_volumes);
+        imported.volumes = import_volumes(*world_, selected.unique_volumes);
         if (selected.particles != DataSelection::none)
         {
             imported.trans_params = import_trans_parameters(selected.particles);
@@ -1256,7 +1238,7 @@ ImportData GeantImporter::operator()(DataSelection const& selected)
  * Return a populated \c ImportVolume vector.
  */
 std::vector<ImportVolume>
-GeantImporter::import_volumes(bool unique_volumes) const
+import_volumes(G4VPhysicalVolume const& world, bool unique_volumes)
 {
     // Note: if the LV has been purged (i.e. by trying to run multiple
     // geometries in the same execution), the instance ID's won't correspond to
@@ -1305,13 +1287,39 @@ GeantImporter::import_volumes(bool unique_volumes) const
                 volume.name = make_gdml_name(lv);
             }
         },
-        *world_);
+        world);
 
     CELER_LOG(debug) << "Loaded " << count << " volumes with "
                      << (unique_volumes ? "uniquified" : "original")
                      << " names";
     return result;
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Create an import patricle.
+ */
+ImportParticle import_particle(G4ParticleDefinition const& p)
+{
+    ImportParticle result;
+    result.name = p.GetParticleName();
+    result.pdg = p.GetPDGEncoding();
+    result.mass = p.GetPDGMass();
+    result.charge = p.GetPDGCharge();
+    result.spin = p.GetPDGSpin();
+    result.lifetime = p.GetPDGLifeTime();
+    result.is_stable = p.GetPDGStable();
+
+    if (!result.is_stable)
+    {
+        double const time_scale = native_value_from_clhep(ImportUnits::time);
+
+        // Convert lifetime of unstable particles to seconds
+        result.lifetime *= time_scale;
+    }
+
+    return result;
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace celeritas

src/celeritas/ext/GeantImporter.hh

@@ -17,6 +17,7 @@
 
 // Geant4 forward declaration
 class G4VPhysicalVolume;  // IWYU pragma: keep
+class G4ParticleDefinition;  // IWYU pragma: keep
 
 namespace celeritas
 {
@@ -96,11 +97,14 @@ class GeantImporter final : public ImporterInterface
     GeantSetup setup_;
     // World physical volume
     G4VPhysicalVolume const* world_{nullptr};
+};
+
+//---------------------------------------------------------------------------//
 
-    //// HELPER FUNCTIONS ////
+std::vector<ImportVolume>
+import_volumes(G4VPhysicalVolume const& world, bool unique_volumes);
 
-    std::vector<ImportVolume> import_volumes(bool unique_volumes) const;
-};
+ImportParticle import_particle(G4ParticleDefinition const& p);
 
 //---------------------------------------------------------------------------//
 // INLINE DEFINITIONS