Remove multiple scattering from Processes

app/demo-geant-integration/NoFieldAlongStepFactory.cc

@@ -10,8 +10,10 @@
 #include <memory>
 #include <type_traits>
 
+#include "celeritas/em/UrbanMscParams.hh"
 #include "celeritas/global/alongstep/AlongStepGeneralLinearAction.hh"
 #include "celeritas/io/ImportData.hh"
+#include "celeritas/phys/ImportedProcessAdapter.hh"
 
 namespace demo_geant
 {
@@ -22,12 +24,12 @@ namespace demo_geant
 auto NoFieldAlongStepFactory::operator()(argument_type input) const
     -> result_type
 {
-    // Create along-step action
     return celeritas::AlongStepGeneralLinearAction::from_params(
         input.action_id,
         *input.material,
         *input.particle,
-        *input.physics,
+        celeritas::UrbanMscParams::from_import(
+            *input.particle, *input.material, *input.imported),
         input.imported->em_params.energy_loss_fluct);
 }
 

app/demo-loop/LDemoIO.cc

@@ -16,6 +16,7 @@
 #include "corecel/io/StringUtils.hh"
 #include "corecel/sys/Device.hh"
 #include "celeritas/Units.hh"
+#include "celeritas/em/UrbanMscParams.hh"
 #include "celeritas/ext/GeantImporter.hh"
 #include "celeritas/ext/GeantPhysicsOptionsIO.json.hh"
 #include "celeritas/ext/GeantSetup.hh"
@@ -272,14 +273,16 @@ TransporterInput load_input(LDemoArgs const& args)
     }
 
     bool eloss = imported_data.em_params.energy_loss_fluct;
+    auto msc = UrbanMscParams::from_import(
+        *params.particle, *params.material, imported_data);
     if (args.mag_field == LDemoArgs::no_field())
     {
         // Create along-step action
         auto along_step = AlongStepGeneralLinearAction::from_params(
             params.action_reg->next_id(),
             *params.material,
             *params.particle,
-            *params.physics,
+            msc,
             eloss);
         params.action_reg->insert(along_step);
     }
@@ -298,8 +301,8 @@ TransporterInput load_input(LDemoArgs const& args)
             f *= units::tesla;
         }
 
-        auto along_step = AlongStepUniformMscAction::from_params(
-            params.action_reg->next_id(), *params.physics, field_params);
+        auto along_step = std::make_shared<AlongStepUniformMscAction>(
+            params.action_reg->next_id(), field_params, msc);
         CELER_ASSERT(along_step->field() != LDemoArgs::no_field());
         params.action_reg->insert(along_step);
     }

example/accel/accel.cc

@@ -39,6 +39,7 @@
 #include <accel/Logger.hh>
 #include <accel/SetupOptions.hh>
 #include <accel/SharedParams.hh>
+#include <celeritas/em/UrbanMscParams.hh>
 #include <celeritas/global/alongstep/AlongStepGeneralLinearAction.hh>
 #include <celeritas/io/ImportData.hh>
 #include <corecel/Macros.hh>
@@ -64,7 +65,8 @@ make_nofield_along_step(celeritas::AlongStepFactoryInput const& input)
         input.action_id,
         *input.material,
         *input.particle,
-        *input.physics,
+        celeritas::UrbanMscParams::from_import(
+            *input.particle, *input.material, *input.imported),
         input.imported->em_params.energy_loss_fluct);
 }
 

src/celeritas/CMakeLists.txt

@@ -29,6 +29,7 @@ list(APPEND SOURCES
   Types.cc
   em/AtomicRelaxationParams.cc
   em/FluctuationParams.cc
+  em/UrbanMscParams.cc
   em/detail/Utils.cc
   em/model/BetheHeitlerModel.cc
   em/model/CombinedBremModel.cc
@@ -40,13 +41,11 @@ list(APPEND SOURCES
   em/model/RayleighModel.cc
   em/model/RelativisticBremModel.cc
   em/model/SeltzerBergerModel.cc
-  em/model/UrbanMscModel.cc
   em/process/BremsstrahlungProcess.cc
   em/process/ComptonProcess.cc
   em/process/EIonizationProcess.cc
   em/process/EPlusAnnihilationProcess.cc
   em/process/GammaConversionProcess.cc
-  em/process/MultipleScatteringProcess.cc
   em/process/PhotoelectricProcess.cc
   em/process/RayleighProcess.cc
   geo/GeoMaterialParams.cc

src/celeritas/Types.hh

@@ -144,6 +144,15 @@ struct NoData
 {
 };
 
+//! Within-kernel state for step limits and action
+struct AlongStepLocalState
+{
+    //! Final step limit (true path length) and action
+    StepLimit step_limit;
+    //! Smooth (straight/curved) step movement without accounting for MSC
+    real_type geo_step{};
+};
+
 //---------------------------------------------------------------------------//
 // HELPER FUNCTIONS (HOST)
 //---------------------------------------------------------------------------//

src/celeritas/em/model/UrbanMscModel.cc → src/celeritas/em/UrbanMscParams.cc

@@ -1,12 +1,13 @@
 //----------------------------------*-C++-*----------------------------------//
-// Copyright 2020-2023 UT-Battelle, LLC, and other Celeritas developers.
+// Copyright 2023 UT-Battelle, LLC, and other Celeritas developers.
 // See the top-level COPYRIGHT file for details.
 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
 //---------------------------------------------------------------------------//
-//! \file celeritas/em/model/UrbanMscModel.cc
+//! \file celeritas/em/UrbanMscParams.cc
 //---------------------------------------------------------------------------//
-#include "UrbanMscModel.hh"
+#include "UrbanMscParams.hh"
 
+#include <algorithm>
 #include <cmath>
 #include <memory>
 #include <utility>
@@ -24,6 +25,7 @@
 #include "celeritas/grid/ValueGridBuilder.hh"
 #include "celeritas/grid/ValueGridInserter.hh"
 #include "celeritas/grid/XsGridData.hh"
+#include "celeritas/io/ImportData.hh"
 #include "celeritas/io/ImportProcess.hh"
 #include "celeritas/mat/MaterialParams.hh"
 #include "celeritas/mat/MaterialView.hh"
@@ -36,22 +38,43 @@ namespace celeritas
 {
 //---------------------------------------------------------------------------//
 /*!
- * Construct from model ID and other necessary data.
+ * Construct if Urban model is present, or else return nullptr.
  */
-UrbanMscModel::UrbanMscModel(ActionId id,
-                             ParticleParams const& particles,
-                             MaterialParams const& materials,
-                             ImportedProcessAdapter const& pdata)
+std::shared_ptr<UrbanMscParams>
+UrbanMscParams::from_import(ParticleParams const& particles,
+                            MaterialParams const& materials,
+                            ImportData const& import)
 {
-    CELER_EXPECT(id);
-    HostValue host_data;
+    auto is_urban = [](ImportMscModel const& imm) {
+        return imm.model == ImportModelClass::urban_msc;
+    };
+    if (!std::any_of(
+            import.msc_models.begin(), import.msc_models.end(), is_urban))
+    {
+        // No Urban MSC present
+        return nullptr;
+    }
+    return std::make_shared<UrbanMscParams>(
+        particles, materials, import.msc_models);
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Construct from cross section data and material properties.
+ */
+UrbanMscParams::UrbanMscParams(ParticleParams const& particles,
+                               MaterialParams const& materials,
+                               VecImportMscModel const& mdata_vec)
+{
+    HostVal<UrbanMscData> host_data;
 
-    host_data.ids.action = id;
     host_data.ids.electron = particles.find(pdg::electron());
     host_data.ids.positron = particles.find(pdg::positron());
     CELER_VALIDATE(host_data.ids.electron && host_data.ids.positron,
-                   << "missing e-/e+ (required for " << this->description()
-                   << ")");
+                   << "missing e-/e+ (required for Urban MSC)");
+
+    // Save electron mass
+    host_data.electron_mass = particles.get(host_data.ids.electron).mass();
 
     // TODO: change IDs to a vector for all particles. This model should apply
     // to muons and charged hadrons as well
@@ -62,16 +85,44 @@ UrbanMscModel::UrbanMscModel(ActionId id,
                               "particles other than electron and positron";
     }
 
-    // Save electron mass
-    host_data.electron_mass = particles.get(host_data.ids.electron).mass();
+    // Filter MSC data by model and particle type
+    std::vector<ImportMscModel const*> urban_data(particles.size(), nullptr);
+    for (ImportMscModel const& imm : mdata_vec)
+    {
+        // Filter out other MSC models
+        if (imm.model != ImportModelClass::urban_msc)
+            continue;
+
+        // Filter out unused particles
+        PDGNumber pdg{imm.particle_pdg};
+        ParticleId pid = pdg ? particles.find(pdg) : ParticleId{};
+        if (!pid)
+            continue;
+
+        // Make sure the model data is unique
+        CELER_VALIDATE(!urban_data[pid.get()],
+                       << "duplicate " << to_cstring(imm.model)
+                       << " physics data for particle "
+                       << particles.id_to_label(pid));
+        urban_data[pid.get()] = &imm;
+    }
+
+    auto get_lambda = [&urban_data, &particles](ParticleId pid) {
+        CELER_ASSERT(pid < urban_data.size());
+        ImportMscModel const* imm = urban_data[pid.unchecked_get()];
+        CELER_VALIDATE(imm,
+                       << "missing Urban MSC physics data for particle "
+                       << particles.id_to_label(pid));
+        return &imm->lambda_table;
+    };
 
     {
         // Particle-dependent data
         Array<ParticleId, 2> const par_ids{
             {host_data.ids.electron, host_data.ids.positron}};
         Array<ImportPhysicsTable const*, 2> const xs_tables{{
-            &pdata.get_lambda(par_ids[0]),
-            &pdata.get_lambda(par_ids[1]),
+            get_lambda(par_ids[0]),
+            get_lambda(par_ids[1]),
         }};
         CELER_ASSERT(xs_tables[0]->x_units == ImportUnits::mev);
         CELER_ASSERT(xs_tables[0]->y_units == ImportUnits::cm_inv);
@@ -97,7 +148,7 @@ UrbanMscModel::UrbanMscModel(ActionId id,
             auto&& mat = materials.get(mat_id);
 
             // Build material-dependent data
-            mdata.push_back(UrbanMscModel::calc_material_data(mat));
+            mdata.push_back(UrbanMscParams::calc_material_data(mat));
 
             // Build particle-dependent data
             const real_type zeff = mat.zeff();
@@ -135,51 +186,6 @@ UrbanMscModel::UrbanMscModel(ActionId id,
     CELER_ENSURE(this->mirror_);
 }
 
-//---------------------------------------------------------------------------//
-/*!
- * Particle types and energy ranges that this model applies to.
- */
-auto UrbanMscModel::applicability() const -> SetApplicability
-{
-    Applicability electron_msc;
-    electron_msc.particle = this->host_ref().ids.electron;
-    electron_msc.lower = zero_quantity();
-    electron_msc.upper = units::MevEnergy{1e+8};
-
-    Applicability positron_msc = electron_msc;
-    positron_msc.particle = this->host_ref().ids.positron;
-
-    return {electron_msc, positron_msc};
-}
-
-//---------------------------------------------------------------------------//
-/*!
- * Get the microscopic cross sections for the given particle and material.
- */
-auto UrbanMscModel::micro_xs(Applicability) const -> MicroXsBuilders
-{
-    // No cross sections for multiple scattering
-    return {};
-}
-
-//---------------------------------------------------------------------------//
-//!@{
-/*!
- * No discrete interaction: it's integrated into along_step.
- */
-void UrbanMscModel::execute(CoreDeviceRef const&) const {}
-void UrbanMscModel::execute(CoreHostRef const&) const {}
-//!@}
-
-//---------------------------------------------------------------------------//
-/*!
- * Get the model ID for this model.
- */
-ActionId UrbanMscModel::action_id() const
-{
-    return this->host_ref().ids.action;
-}
-
 //---------------------------------------------------------------------------//
 /*!
  * Build UrbanMsc data per material.
@@ -188,11 +194,11 @@ ActionId UrbanMscModel::action_id() const
  * documented in section 8.1.5 of the Geant4 10.7 Physics Reference Manual.
  */
 UrbanMscMaterialData
-UrbanMscModel::calc_material_data(MaterialView const& material_view)
+UrbanMscParams::calc_material_data(MaterialView const& material_view)
 {
     using PolyQuad = PolyEvaluator<double, 2>;
 
-    MaterialData data;
+    UrbanMscMaterialData data;
 
     double zeff = material_view.zeff();
 
@@ -219,6 +225,5 @@ UrbanMscModel::calc_material_data(MaterialView const& material_view)
 
     return data;
 }
-
 //---------------------------------------------------------------------------//
 }  // namespace celeritas