Import a few more EM parameters from Geant4

app/celer-dump-data.cc

@@ -408,20 +408,23 @@ void print_em_params(ImportEmParameters const& em_params)
     // NOTE: boolalpha doesn't work with setw, see
     // https://timsong-cpp.github.io/lwg-issues/2703
 #define PEP_STREAM_PARAM(NAME) \
-    "| " << setw(18) << #NAME << " | " << setw(7) << em_params.NAME << " |\n"
+    "| " << setw(22) << #NAME << " | " << setw(7) << em_params.NAME << " |\n"
 #define PEP_STREAM_BOOL(NAME)                     \
-    "| " << setw(18) << #NAME << " | " << setw(7) \
+    "| " << setw(22) << #NAME << " | " << setw(7) \
          << (em_params.NAME ? "true" : "false") << " |\n"
     cout << R"gfm(
 # EM parameters
 
-| EM parameter       | Value   |
-| ------------------ | ------- |
+| EM parameter           | Value   |
+| ---------------------- | ------- |
 )gfm";
     cout << PEP_STREAM_BOOL(energy_loss_fluct) << PEP_STREAM_BOOL(lpm)
          << PEP_STREAM_BOOL(integral_approach)
-         << PEP_STREAM_PARAM(linear_loss_limit) << PEP_STREAM_BOOL(auger)
-         << endl;
+         << PEP_STREAM_PARAM(linear_loss_limit)
+         << PEP_STREAM_PARAM(lowest_electron_energy) << PEP_STREAM_BOOL(auger)
+         << PEP_STREAM_PARAM(msc_range_factor)
+         << PEP_STREAM_PARAM(msc_safety_factor)
+         << PEP_STREAM_PARAM(msc_lambda_limit) << endl;
 #undef PEP_STREAM_PARAM
 #undef PEP_STREAM_BOOL
 }

app/demo-loop/LDemoIO.cc

@@ -272,6 +272,8 @@ TransporterInput load_input(LDemoArgs const& args)
         input.options.fixed_step_limiter = args.step_limiter;
         input.options.secondary_stack_factor = args.secondary_stack_factor;
         input.options.linear_loss_limit = imported.em_params.linear_loss_limit;
+        input.options.lowest_electron_energy = PhysicsParamsOptions::Energy{
+            imported.em_params.lowest_electron_energy};
 
         input.processes = [&params, &args, &imported] {
             std::vector<std::shared_ptr<Process const>> result;

src/accel/SharedParams.cc

@@ -325,6 +325,8 @@ void SharedParams::initialize_core(SetupOptions const& options)
         input.action_registry = params.action_reg.get();
 
         input.options.linear_loss_limit = imported->em_params.linear_loss_limit;
+        input.options.lowest_electron_energy = PhysicsParamsOptions::Energy{
+            imported->em_params.lowest_electron_energy};
         input.options.secondary_stack_factor = options.secondary_stack_factor;
 
         return std::make_shared<PhysicsParams>(std::move(input));

src/celeritas/em/UrbanMscParams.cc

@@ -54,8 +54,14 @@ UrbanMscParams::from_import(ParticleParams const& particles,
         // No Urban MSC present
         return nullptr;
     }
+
+    Options opts;
+    opts.lambda_limit = import.em_params.msc_lambda_limit;
+    opts.safety_fact = import.em_params.msc_safety_factor;
+    opts.range_fact = import.em_params.msc_range_factor;
+
     return std::make_shared<UrbanMscParams>(
-        particles, materials, import.msc_models);
+        particles, materials, import.msc_models, opts);
 }
 
 //---------------------------------------------------------------------------//
@@ -64,7 +70,8 @@ UrbanMscParams::from_import(ParticleParams const& particles,
  */
 UrbanMscParams::UrbanMscParams(ParticleParams const& particles,
                                MaterialParams const& materials,
-                               VecImportMscModel const& mdata_vec)
+                               VecImportMscModel const& mdata_vec,
+                               Options options)
 {
     HostVal<UrbanMscData> host_data;
 
@@ -85,6 +92,20 @@ UrbanMscParams::UrbanMscParams(ParticleParams const& particles,
                               "particles other than electron and positron";
     }
 
+    // Save parameters
+    CELER_VALIDATE(options.lambda_limit > 0,
+                   << "invalid lambda_limit=" << options.lambda_limit
+                   << " (should be positive)");
+    CELER_VALIDATE(options.safety_fact >= 0.1,
+                   << "invalid safety_fact=" << options.safety_fact
+                   << " (should be >= 0.1)");
+    CELER_VALIDATE(options.range_fact > 0 && options.range_fact < 1,
+                   << "invalid range_fact=" << options.range_fact
+                   << " (should be within 0 < limit < 1)");
+    host_data.params.lambda_limit = options.lambda_limit;
+    host_data.params.range_fact = options.range_fact;
+    host_data.params.safety_fact = options.safety_fact;
+
     // Filter MSC data by model and particle type
     std::vector<ImportMscModel const*> urban_data(particles.size(), nullptr);
     for (ImportMscModel const& imm : mdata_vec)

src/celeritas/em/UrbanMscParams.hh

@@ -42,6 +42,14 @@ class UrbanMscParams
     using VecImportMscModel = std::vector<ImportMscModel>;
     //!@}
 
+    //! MSC configuration options
+    struct Options
+    {
+        real_type lambda_limit{1 * units::millimeter};  //!< Lambda limit
+        real_type range_fact{0.04};  //!< Range factor for e-/e+
+        real_type safety_fact{0.6};  //!< Safety factor
+    };
+
   public:
     // Construct if MSC process data is present, else return nullptr
     static std::shared_ptr<UrbanMscParams>
@@ -52,7 +60,8 @@ class UrbanMscParams
     // Construct from process data
     inline UrbanMscParams(ParticleParams const& particles,
                           MaterialParams const& materials,
-                          VecImportMscModel const& mdata);
+                          VecImportMscModel const& mdata,
+                          Options options);
 
     // TODO: possible "applicability" interface used for constructing
     // along-step kernels?

src/celeritas/ext/GeantImporter.cc

@@ -600,7 +600,11 @@ ImportEmParameters import_em_parameters()
     import.lpm = g4.LPM();
     import.integral_approach = g4.Integral();
     import.linear_loss_limit = g4.LinearLossLimit();
+    import.lowest_electron_energy = g4.LowestElectronEnergy() / MeV;
     import.auger = g4.Auger();
+    import.msc_range_factor = g4.MscRangeFactor();
+    import.msc_safety_factor = g4.MscSafetyFactor();
+    import.msc_lambda_limit = g4.MscLambdaLimit() / cm;
 
     CELER_ENSURE(import);
     return import;

src/celeritas/ext/GeantPhysicsOptions.hh

@@ -61,6 +61,10 @@ enum class RelaxationSelection
  * - \c min_energy: lowest energy of any EM physics process
  * - \c max_energy: highest energy of any EM physics process
  * - \c linear_loss_limit: see \c PhysicsParamsOptions::linear_loss_limit
+ * - \c lowest_electron_energy: lowest e-/e+ kinetic energy
+ * - \c msc_range_factor: e-/e+ range factor for MSC models
+ * - \c msc_safety_factor: safety factor for MSC models
+ * - \c msc_lambda_limit: lambda limit for MSC models [cm]
  * - \c verbose: print detailed Geant4 output
  */
 struct GeantPhysicsOptions
@@ -80,6 +84,10 @@ struct GeantPhysicsOptions
     units::MevEnergy min_energy{0.1 * 1e-3};  // 0.1 keV
     units::MevEnergy max_energy{100 * 1e6};  // 100 TeV
     real_type linear_loss_limit{0.01};
+    units::MevEnergy lowest_electron_energy{0.001};  // 1 keV
+    real_type msc_range_factor{0.04};
+    real_type msc_safety_factor{0.6};
+    real_type msc_lambda_limit{0.1};  // 1 mm
 
     bool verbose{false};
 };
@@ -98,6 +106,10 @@ operator==(GeantPhysicsOptions const& a, GeantPhysicsOptions const& b)
            && a.em_bins_per_decade == b.em_bins_per_decade
            && a.min_energy == b.min_energy && a.max_energy == b.max_energy
            && a.linear_loss_limit == b.linear_loss_limit
+           && a.lowest_electron_energy == b.lowest_electron_energy
+           && a.msc_range_factor == b.msc_range_factor
+           && a.msc_safety_factor == b.msc_safety_factor
+           && a.msc_lambda_limit == b.msc_lambda_limit
            && a.verbose == b.verbose;
 }
 

src/celeritas/ext/GeantPhysicsOptionsIO.json.cc

@@ -136,6 +136,10 @@ void from_json(nlohmann::json const& j, GeantPhysicsOptions& options)
     GPO_LOAD_OPTION(min_energy);
     GPO_LOAD_OPTION(max_energy);
     GPO_LOAD_OPTION(linear_loss_limit);
+    GPO_LOAD_OPTION(lowest_electron_energy);
+    GPO_LOAD_OPTION(msc_range_factor);
+    GPO_LOAD_OPTION(msc_safety_factor);
+    GPO_LOAD_OPTION(msc_lambda_limit);
     GPO_LOAD_OPTION(verbose);
 #undef GPO_LOAD_OPTION
 }
@@ -161,6 +165,10 @@ void to_json(nlohmann::json& j, GeantPhysicsOptions const& options)
     GPO_SAVE_OPTION(min_energy);
     GPO_SAVE_OPTION(max_energy);
     GPO_SAVE_OPTION(linear_loss_limit);
+    GPO_SAVE_OPTION(lowest_electron_energy);
+    GPO_SAVE_OPTION(msc_range_factor);
+    GPO_SAVE_OPTION(msc_safety_factor);
+    GPO_SAVE_OPTION(msc_lambda_limit);
     GPO_SAVE_OPTION(verbose);
 #undef GPO_SAVE_OPTION
 }

src/celeritas/ext/detail/GeantPhysicsList.cc

@@ -67,6 +67,12 @@ GeantPhysicsList::GeantPhysicsList(Options const& options) : options_(options)
     em_parameters.SetAuger(options.relaxation == RelaxationSelection::all);
     em_parameters.SetIntegral(options.integral_approach);
     em_parameters.SetLinearLossLimit(options.linear_loss_limit);
+    em_parameters.SetMscRangeFactor(options.msc_range_factor);
+    em_parameters.SetMscSafetyFactor(options.msc_safety_factor);
+    em_parameters.SetMscLambdaLimit(options.msc_lambda_limit * CLHEP::cm);
+    em_parameters.SetLowestElectronEnergy(
+        value_as<units::MevEnergy>(options.lowest_electron_energy)
+        * CLHEP::MeV);
 
     int verb = options_.verbose ? 1 : 0;
     this->SetVerboseLevel(verb);

src/celeritas/global/alongstep/AlongStep.hh

@@ -156,9 +156,10 @@ inline CELER_FUNCTION void along_step(MH&& msc,
             particle.subtract_energy(deposited);
         }
         // Energy loss helper *must* apply the tracking cutoff
-        CELER_ASSERT(particle.energy()
-                         >= track.make_physics_view().scalars().eloss_calc_limit
-                     || !apply_cut || particle.is_stopped());
+        CELER_ASSERT(
+            particle.energy()
+                >= track.make_physics_view().scalars().lowest_electron_energy
+            || !apply_cut || particle.is_stopped());
     }
 
     if (particle.is_stopped())

src/celeritas/global/alongstep/detail/FluctELoss.hh

@@ -102,7 +102,7 @@ CELER_FUNCTION auto FluctELoss::calc_eloss(CoreTrackView const& track,
     auto particle = track.make_particle_view();
     auto phys = track.make_physics_view();
 
-    if (apply_cut && particle.energy() < phys.scalars().eloss_calc_limit)
+    if (apply_cut && particle.energy() < phys.scalars().lowest_electron_energy)
     {
         // Deposit all energy immediately when we start below the tracking cut
         return particle.energy();
@@ -158,7 +158,7 @@ CELER_FUNCTION auto FluctELoss::calc_eloss(CoreTrackView const& track,
     }
 
     if (apply_cut
-        && (particle.energy() - eloss <= phys.scalars().eloss_calc_limit))
+        && (particle.energy() - eloss <= phys.scalars().lowest_electron_energy))
     {
         // Deposit all energy when we end below the tracking cut
         return particle.energy();

src/celeritas/global/alongstep/detail/MeanELoss.hh

@@ -66,7 +66,7 @@ CELER_FUNCTION auto MeanELoss::calc_eloss(CoreTrackView const& track,
     auto particle = track.make_particle_view();
     auto phys = track.make_physics_view();
 
-    if (apply_cut && particle.energy() < phys.scalars().eloss_calc_limit)
+    if (apply_cut && particle.energy() < phys.scalars().lowest_electron_energy)
     {
         // Deposit all energy when we start below the tracking cut
         return particle.energy();
@@ -76,7 +76,7 @@ CELER_FUNCTION auto MeanELoss::calc_eloss(CoreTrackView const& track,
     Energy eloss = calc_mean_energy_loss(particle, phys, step);
 
     if (apply_cut
-        && (particle.energy() - eloss <= phys.scalars().eloss_calc_limit))
+        && (particle.energy() - eloss <= phys.scalars().lowest_electron_energy))
     {
         // Deposit all energy when we end below the tracking cut
         return particle.energy();

src/celeritas/io/ImportParameters.hh

@@ -16,6 +16,10 @@ namespace celeritas
  * Common electromagnetic physics parameters (see G4EmParameters.hh).
  *
  * \note Geant4 v11 removed the Spline() option from G4EmParameters.hh.
+ * \note The Geant4 MSC models use the values in \c G4EmParameters as the
+ * defaults; however, the MSC parameters can also be set directly using the
+ * model setter methods (there is no way to retrieve the values from the model
+ * in that case).
  */
 struct ImportEmParameters
 {
@@ -27,11 +31,24 @@ struct ImportEmParameters
     bool integral_approach{true};
     //! Slowing down threshold for linearity assumption
     double linear_loss_limit{0.01};
+    //! Lowest e-/e+ kinetic energy [MeV]
+    double lowest_electron_energy{0.001};
     //! Whether auger emission should be enabled (valid only for relaxation)
     bool auger{false};
+    //! MSC range factor for e-/e+
+    double msc_range_factor{0.04};
+    //! MSC safety factor
+    double msc_safety_factor{0.6};
+    //! MSC lambda limit [cm]
+    double msc_lambda_limit{0.1};
 
     //! Whether parameters are assigned and valid
-    explicit operator bool() const { return linear_loss_limit > 0; }
+    explicit operator bool() const
+    {
+        return linear_loss_limit > 0 && lowest_electron_energy > 0
+               && msc_range_factor > 0 && msc_range_factor < 1
+               && msc_safety_factor >= 0.1 && msc_lambda_limit > 0;
+    }
 };
 
 //---------------------------------------------------------------------------//

src/celeritas/phys/PhysicsData.hh

@@ -222,7 +222,7 @@ struct PhysicsParamsScalars
     real_type min_range{};  //!< rho [cm]
     real_type max_step_over_range{};  //!< alpha [unitless]
     real_type min_eprime_over_e{};  //!< xi [unitless]
-    Energy eloss_calc_limit{};  //!< Lowest energy for eloss calculation
+    Energy lowest_electron_energy{};  //!< Lowest e-/e+ kinetic energy
     real_type linear_loss_limit{};  //!< For scaled range calculation
     real_type fixed_step_limiter{};  //!< Global charged step size limit [cm]
 
@@ -237,7 +237,8 @@ struct PhysicsParamsScalars
     {
         return max_particle_processes > 0 && model_to_action >= 4
                && num_models > 0 && min_range > 0 && max_step_over_range > 0
-               && min_eprime_over_e > 0 && eloss_calc_limit > zero_quantity()
+               && min_eprime_over_e > 0
+               && lowest_electron_energy > zero_quantity()
                && linear_loss_limit > 0 && secondary_stack_factor > 0
                && ((fixed_step_limiter > 0)
                    == static_cast<bool>(fixed_step_action));

src/celeritas/phys/PhysicsParams.cc

@@ -220,9 +220,10 @@ void PhysicsParams::build_options(Options const& opts, HostValue* data) const
     CELER_VALIDATE(opts.min_range > 0,
                    << "invalid min_range=" << opts.min_range
                    << " (should be positive)");
-    CELER_VALIDATE(opts.eloss_calc_limit.value() > 0,
-                   << "invalid eloss_calc_limit="
-                   << opts.eloss_calc_limit.value() << " (should be positive)");
+    CELER_VALIDATE(opts.lowest_electron_energy.value() > 0,
+                   << "invalid lowest_electron_energy="
+                   << opts.lowest_electron_energy.value()
+                   << " (should be positive)");
     CELER_VALIDATE(opts.linear_loss_limit >= 0 && opts.linear_loss_limit <= 1,
                    << "invalid linear_loss_limit=" << opts.linear_loss_limit
                    << " (should be within 0 <= limit <= 1)");
@@ -232,7 +233,7 @@ void PhysicsParams::build_options(Options const& opts, HostValue* data) const
     data->scalars.min_range = opts.min_range;
     data->scalars.max_step_over_range = opts.max_step_over_range;
     data->scalars.min_eprime_over_e = opts.min_eprime_over_e;
-    data->scalars.eloss_calc_limit = opts.eloss_calc_limit;
+    data->scalars.lowest_electron_energy = opts.lowest_electron_energy;
     data->scalars.linear_loss_limit = opts.linear_loss_limit;
     data->scalars.secondary_stack_factor = opts.secondary_stack_factor;
 }

src/celeritas/phys/PhysicsParams.hh

@@ -50,8 +50,7 @@ class ParticleParams;
  * - \c linear_loss_limit: if the mean energy loss along a step is greater than
  *   this fractional value of the pre-step kinetic energy, recalculate the
  *   energy loss.
- * - \c eloss_calc_limit: kill charged particles below this energy at the end
- *   of a step.
+ * - \c lowest_electron_energy: lowest kinetic energy for electrons/positrons
  * - \c secondary_stack_factor: the number of secondary slots per track slot
  *   allocated.
  * - \c disable_integral_xs: for particles with energy loss processes, the
@@ -62,7 +61,8 @@ class ParticleParams;
  *   processes.
  *
  * NOTE: min_range/max_step_over_range are not accessible through Geant4, and
- * they can also be set to be different for electrons, mu/hadrons, and ions.
+ * they can also be set to be different for electrons, mu/hadrons, and ions
+ * (they are set in Geant4 with \c G4EmParameters::SetStepFunction()).
  */
 struct PhysicsParamsOptions
 {
@@ -79,7 +79,7 @@ struct PhysicsParamsOptions
     //! \name Energy loss
     real_type min_eprime_over_e = 0.8;
     real_type linear_loss_limit = 0.01;
-    Energy eloss_calc_limit = Energy{0.001};
+    Energy lowest_electron_energy = Energy{0.001};
     //!@}
 
     real_type secondary_stack_factor = 3;

src/celeritas/phys/PhysicsParamsOutput.cc

@@ -84,7 +84,7 @@ void PhysicsParamsOutput::output(JsonPimpl* j) const
         PPO_SAVE_OPTION(min_range);
         PPO_SAVE_OPTION(max_step_over_range);
         PPO_SAVE_OPTION(min_eprime_over_e);
-        PPO_SAVE_OPTION(eloss_calc_limit);
+        PPO_SAVE_OPTION(lowest_electron_energy);
         PPO_SAVE_OPTION(linear_loss_limit);
         PPO_SAVE_OPTION(fixed_step_limiter);
 #    undef PPO_SAVE_OPTION

test/celeritas/data/four-steel-slabs.geant.json

@@ -5,6 +5,10 @@
 "em_bins_per_decade": 7,
 "integral_approach": true,
 "linear_loss_limit": 0.01,
+"lowest_electron_energy": 0.001,
+"msc_range_factor": 0.04,
+"msc_safety_factor": 0.6,
+"msc_lambda_limit": 0.1,
 "lpm": true,
 "max_energy": [
  1000000.0,