diff --git a/app/celer-sim/Runner.cc b/app/celer-sim/Runner.cc index 8e0e37ffab..2d480af887 100644 --- a/app/celer-sim/Runner.cc +++ b/app/celer-sim/Runner.cc @@ -30,6 +30,7 @@ #include "celeritas/Types.hh" #include "celeritas/Units.hh" #include "celeritas/em/params/UrbanMscParams.hh" +#include "celeritas/em/params/WentzelOKVIParams.hh" #include "celeritas/ext/GeantImporter.hh" #include "celeritas/ext/GeantSetup.hh" #include "celeritas/ext/RootFileManager.hh" @@ -308,6 +309,9 @@ void Runner::build_core_params(RunnerInput const& inp, params.cutoff = CutoffParams::from_import( imported, params.particle, params.material); + // Construct shared data for Coulomb scattering + params.wentzel = WentzelOKVIParams::from_import(imported, params.material); + // Load physics: create individual processes with make_shared params.physics = [¶ms, &inp, &imported] { PhysicsParams::Input input; diff --git a/src/accel/SharedParams.cc b/src/accel/SharedParams.cc index 0a882de325..ea6483fd4f 100644 --- a/src/accel/SharedParams.cc +++ b/src/accel/SharedParams.cc @@ -36,6 +36,7 @@ #include "geocel/GeantUtils.hh" #include "geocel/g4/GeantGeoParams.hh" #include "celeritas/Types.hh" +#include "celeritas/em/params/WentzelOKVIParams.hh" #include "celeritas/ext/GeantImporter.hh" #include "celeritas/ext/RootExporter.hh" #include "celeritas/geo/GeoMaterialParams.hh" @@ -58,6 +59,7 @@ #include "AlongStepFactory.hh" #include "SetupOptions.hh" + #include "detail/HitManager.hh" #include "detail/OffloadWriter.hh" @@ -530,6 +532,9 @@ void SharedParams::initialize_core(SetupOptions const& options) params.cutoff = CutoffParams::from_import( *imported, params.particle, params.material); + // Construct shared data for Coulomb scattering + params.wentzel = WentzelOKVIParams::from_import(*imported, params.material); + // Load physics: create individual processes with make_shared params.physics = [¶ms, &options, &imported] { PhysicsParams::Input input; diff --git a/src/celeritas/CMakeLists.txt b/src/celeritas/CMakeLists.txt index 87710fdf05..04084ddd52 100644 --- a/src/celeritas/CMakeLists.txt +++ b/src/celeritas/CMakeLists.txt @@ -18,6 +18,7 @@ list(APPEND SOURCES em/params/AtomicRelaxationParams.cc em/params/FluctuationParams.cc em/params/UrbanMscParams.cc + em/params/WentzelOKVIParams.cc em/params/WentzelVIMscParams.cc em/params/detail/MscParamsHelper.cc em/process/BremsstrahlungProcess.cc diff --git a/src/celeritas/Types.cc b/src/celeritas/Types.cc index a1616dcb39..398256a802 100644 --- a/src/celeritas/Types.cc +++ b/src/celeritas/Types.cc @@ -93,6 +93,21 @@ char const* to_cstring(MscStepLimitAlgorithm value) return to_cstring_impl(value); } +//---------------------------------------------------------------------------// +/*! + * Get a string corresponding to the nuclear form factor model. + */ +char const* to_cstring(NuclearFormFactorType value) +{ + static EnumStringMapper const to_cstring_impl{ + "none", + "flat", + "exponential", + "gaussian", + }; + return to_cstring_impl(value); +} + //---------------------------------------------------------------------------// /*! * Checks that the TrackOrder will sort tracks by actions applied at the given diff --git a/src/celeritas/Types.hh b/src/celeritas/Types.hh index 9b154de233..8fb6377dea 100644 --- a/src/celeritas/Types.hh +++ b/src/celeritas/Types.hh @@ -160,6 +160,17 @@ enum class MscStepLimitAlgorithm size_, }; +//---------------------------------------------------------------------------// +//! Nuclear form factor model for Coulomb scattering +enum class NuclearFormFactorType +{ + none, + flat, + exponential, + gaussian, + size_ +}; + //---------------------------------------------------------------------------// // HELPER STRUCTS //---------------------------------------------------------------------------// @@ -196,6 +207,9 @@ char const* to_cstring(TrackOrder); // Get a string corresponding to the MSC step limit algorithm char const* to_cstring(MscStepLimitAlgorithm value); +// Get a string corresponding to the nuclear form factor model +char const* to_cstring(NuclearFormFactorType value); + // Checks that the TrackOrder will sort tracks by actions applied at the given // ActionOrder bool is_action_sorted(ActionOrder action, TrackOrder track); diff --git a/src/celeritas/em/data/MscData.hh b/src/celeritas/em/data/CommonCoulombData.hh similarity index 93% rename from src/celeritas/em/data/MscData.hh rename to src/celeritas/em/data/CommonCoulombData.hh index bab60d6930..e8ae40612f 100644 --- a/src/celeritas/em/data/MscData.hh +++ b/src/celeritas/em/data/CommonCoulombData.hh @@ -3,7 +3,7 @@ // See the top-level COPYRIGHT file for details. // SPDX-License-Identifier: (Apache-2.0 OR MIT) //---------------------------------------------------------------------------// -//! \file celeritas/em/data/MscData.hh +//! \file celeritas/em/data/CommonCoulombData.hh //---------------------------------------------------------------------------// #pragma once @@ -19,7 +19,7 @@ namespace celeritas * * TODO these will probably be changed to a map over all particle IDs. */ -struct MscIds +struct CoulombIds { ParticleId electron; ParticleId positron; diff --git a/src/celeritas/em/data/CoulombScatteringData.hh b/src/celeritas/em/data/CoulombScatteringData.hh index 5c068beda7..220fea4de1 100644 --- a/src/celeritas/em/data/CoulombScatteringData.hh +++ b/src/celeritas/em/data/CoulombScatteringData.hh @@ -9,120 +9,31 @@ #include "corecel/Macros.hh" #include "corecel/Types.hh" -#include "corecel/data/Collection.hh" #include "celeritas/Quantities.hh" #include "celeritas/Types.hh" -namespace celeritas -{ -//---------------------------------------------------------------------------// -/*! - * Particle and action ids used by CoulombScatteringModel. - */ -struct CoulombScatteringIds -{ - ActionId action; - ParticleId electron; - ParticleId positron; - ParticleId proton; - - explicit CELER_FUNCTION operator bool() const - { - // TODO: enable when protons are supported - return action && electron && positron /* && proton */; - } -}; - -//---------------------------------------------------------------------------// -/*! - * Per-element data used by the CoulombScatteringModel. - * - * The matrix of coefficients used to approximate the ratio of the Mott to - * Rutherford cross sections was developed in - * T. Lijian, H. Quing and L. Zhengming, Radiat. Phys. Chem. 45 (1995), - * 235-245 - * and - * M. J. Boschini et al. arXiv:1111.4042 - */ -struct CoulombScatteringElementData -{ - //!@{ - //! \name Dimensions for Mott coefficient matrices - static constexpr size_type num_mott_beta_bins = 6; - static constexpr size_type num_mott_theta_bins = 5; - static constexpr size_type num_mott_elements = 92; - //!@} - - using BetaArray = Array; - using ThetaArray = Array; - using MottCoeffMatrix = Array; +#include "CommonCoulombData.hh" - //! Matrix of Mott coefficients [theta][beta] - MottCoeffMatrix mott_coeff; -}; - -//---------------------------------------------------------------------------// -/*! - * Supported models of nuclear form factors. - */ -enum class NuclearFormFactorType +namespace celeritas { - none, - flat, - exponential, - gaussian -}; - //---------------------------------------------------------------------------// /*! * Constant shared data used by the CoulombScatteringModel. */ -template struct CoulombScatteringData { - template - using ElementItems = celeritas::Collection; - template - using IsotopeItems = celeritas::Collection; - - // Ids - CoulombScatteringIds ids; - - //! Constant prefactor for the squared momentum transfer [(MeV/c)^-2] - IsotopeItems nuclear_form_prefactor; - - // Per element form factors - ElementItems elem_data; + // Particle IDs + CoulombIds ids; - // User-defined factor for the screening coefficient - real_type screening_factor; + // Action ID + ActionId action; - // Model for the form factor to use - NuclearFormFactorType form_factor_type; + //! Cosine of the maximum scattering polar angle + static CELER_CONSTEXPR_FUNCTION real_type cos_thetamax() { return -1; } // Check if the data is initialized - explicit CELER_FUNCTION operator bool() const - { - return ids && !nuclear_form_prefactor.empty() && !elem_data.empty(); - } - - // Copy initialize from an existing CoulombScatteringData - template - CoulombScatteringData& operator=(CoulombScatteringData const& other) - { - CELER_EXPECT(other); - ids = other.ids; - nuclear_form_prefactor = other.nuclear_form_prefactor; - elem_data = other.elem_data; - form_factor_type = other.form_factor_type; - screening_factor = other.screening_factor; - return *this; - } + explicit CELER_FUNCTION operator bool() const { return ids && action; } }; -using CoulombScatteringDeviceRef = DeviceCRef; -using CoulombScatteringHostRef = HostCRef; -using CoulombScatteringRef = NativeCRef; - //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/em/data/UrbanMscData.hh b/src/celeritas/em/data/UrbanMscData.hh index 1c28a90aaf..1e6873fb38 100644 --- a/src/celeritas/em/data/UrbanMscData.hh +++ b/src/celeritas/em/data/UrbanMscData.hh @@ -14,7 +14,7 @@ #include "celeritas/Types.hh" #include "celeritas/grid/XsGridData.hh" -#include "MscData.hh" +#include "CommonCoulombData.hh" namespace celeritas { @@ -137,7 +137,7 @@ struct UrbanMscData //// DATA //// //! Particle IDs - MscIds ids; + CoulombIds ids; //! Mass of of electron in MeV units::MevMass electron_mass; //! User-assignable options diff --git a/src/celeritas/em/data/WentzelOKVIData.hh b/src/celeritas/em/data/WentzelOKVIData.hh new file mode 100644 index 0000000000..3a87c88629 --- /dev/null +++ b/src/celeritas/em/data/WentzelOKVIData.hh @@ -0,0 +1,123 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 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/data/WentzelOKVIData.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Macros.hh" +#include "corecel/Types.hh" +#include "corecel/data/Collection.hh" +#include "celeritas/Constants.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/Types.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Parameters used in both single Coulomb scattering and Wentzel VI MSC models. + * + * When the single Coulomb scattering and Wentzel VI MSC models are used + * together, the MSC model is used to sample scatterings with angles below the + * polar angle limit, and the single scattering model is used for angles above + * the limit. + */ +struct CoulombParameters +{ + //! Whether to use combined single and multiple scattering + bool is_combined{true}; + //! Polar angle limit between single and multiple scattering + real_type costheta_limit{-1}; + //! Factor for the screening coefficient + real_type screening_factor{1}; + //! Factor used to calculate the maximum scattering angle off of a nucleus + real_type a_sq_factor{0.5 + * ipow<2>(constants::hbar_planck * constants::c_light + * units::femtometer)}; + // Model for the form factor to use + NuclearFormFactorType form_factor_type{NuclearFormFactorType::exponential}; + + explicit CELER_FUNCTION operator bool() const + { + return costheta_limit >= -1 && costheta_limit <= 1 + && screening_factor > 0 && a_sq_factor >= 0 + && form_factor_type != NuclearFormFactorType::size_; + } +}; + +//---------------------------------------------------------------------------// +/*! + * Per-element data used by the Coulomb scattering and Wentzel VI models. + * + * The matrix of coefficients used to approximate the ratio of the Mott to + * Rutherford cross sections was developed in T. Lijian, H. Quing and L. + * Zhengming, Radiat. Phys. Chem. 45 (1995), 235-245 and M. J. Boschini et al. + * arXiv:1111.4042 + */ +struct MottElementData +{ + //!@{ + //! \name Dimensions for Mott coefficient matrices + static constexpr size_type num_mott_beta_bins = 6; + static constexpr size_type num_mott_theta_bins = 5; + static constexpr size_type num_mott_elements = 92; + //!@} + + using BetaArray = Array; + using ThetaArray = Array; + using MottCoeffMatrix = Array; + + //! Matrix of Mott coefficients [theta][beta] + MottCoeffMatrix mott_coeff; +}; + +//---------------------------------------------------------------------------// +/*! + * Constant shared data used by the Coulomb scattering and Wentzel VI models. + */ +template +struct WentzelOKVIData +{ + template + using ElementItems = celeritas::Collection; + template + using IsotopeItems = celeritas::Collection; + template + using MaterialItems = Collection; + + // User-assignable parameters + CoulombParameters params; + + // Constant prefactor for the squared momentum transfer [(MeV/c)^-2] + IsotopeItems nuclear_form_prefactor; + + // Per element form factors + ElementItems elem_data; + + // Inverse effective A^2/3 [1/mass^2/3] + MaterialItems inv_mass_cbrt_sq; + + // Check if the data is initialized + explicit CELER_FUNCTION operator bool() const + { + return params && !elem_data.empty() + && params.is_combined == !inv_mass_cbrt_sq.empty(); + } + + template + WentzelOKVIData& operator=(WentzelOKVIData const& other) + { + CELER_EXPECT(other); + params = other.params; + nuclear_form_prefactor = other.nuclear_form_prefactor; + elem_data = other.elem_data; + inv_mass_cbrt_sq = other.inv_mass_cbrt_sq; + return *this; + } +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/data/WentzelVIMscData.hh b/src/celeritas/em/data/WentzelVIMscData.hh index 0424b8559b..24a54b5a9f 100644 --- a/src/celeritas/em/data/WentzelVIMscData.hh +++ b/src/celeritas/em/data/WentzelVIMscData.hh @@ -14,7 +14,7 @@ #include "celeritas/Types.hh" #include "celeritas/grid/XsGridData.hh" -#include "MscData.hh" +#include "CommonCoulombData.hh" namespace celeritas { @@ -26,7 +26,7 @@ struct WentzelVIMscParameters { using Energy = units::MevEnergy; - real_type single_scattering_fact{1.25}; //!< single scattering factor + real_type single_scattering_factor{1.25}; Energy low_energy_limit{0}; Energy high_energy_limit{0}; @@ -52,7 +52,7 @@ struct WentzelVIMscData //// DATA //// //! Particle IDs - MscIds ids; + CoulombIds ids; //! Mass of of electron in MeV units::MevMass electron_mass; //! User-assignable options diff --git a/src/celeritas/em/distribution/WentzelDistribution.hh b/src/celeritas/em/distribution/WentzelDistribution.hh index fa9dcda109..0c95d80d15 100644 --- a/src/celeritas/em/distribution/WentzelDistribution.hh +++ b/src/celeritas/em/distribution/WentzelDistribution.hh @@ -12,7 +12,6 @@ #include "corecel/Macros.hh" #include "corecel/Types.hh" #include "corecel/math/Algorithms.hh" -#include "celeritas/em/data/CoulombScatteringData.hh" #include "celeritas/em/interactor/detail/PhysicsConstants.hh" #include "celeritas/em/xs/MottRatioCalculator.hh" #include "celeritas/em/xs/WentzelHelper.hh" @@ -53,11 +52,13 @@ class WentzelDistribution public: // Construct with state and model data inline CELER_FUNCTION - WentzelDistribution(ParticleTrackView const& particle, + WentzelDistribution(NativeCRef const& wentzel, + WentzelHelper const& helper, + ParticleTrackView const& particle, IsotopeView const& target, - CoulombScatteringElementData const& element_data, - Energy cutoff, - CoulombScatteringRef const& data); + ElementId el_id, + real_type cos_thetamin, + real_type cos_thetamax); // Sample the polar scattering angle template @@ -66,8 +67,11 @@ class WentzelDistribution private: //// DATA //// - // Shared model data - CoulombScatteringRef const& data_; + // Shared Coulomb scattering data + NativeCRef const& wentzel_; + + // Helper for calculating xs ratio and other quantities + WentzelHelper const& helper_; // Incident particle ParticleTrackView const& particle_; @@ -75,11 +79,16 @@ class WentzelDistribution // Target isotope IsotopeView const& target_; - // Mott coefficients for the target element - CoulombScatteringElementData const& element_data_; + // Target element + ElementId el_id_; - // Helper for calculating xs ratio and other quantities - WentzelHelper const helper_; + // Cosine of the minimum scattering angle + real_type cos_thetamin_; + + // Cosine of the maximum scattering angle + real_type cos_thetamax_; + + //// HELPER FUNCTIONS //// // Calculates the form factor from the scattered polar angle inline CELER_FUNCTION real_type calculate_form_factor(real_type cos_t) const; @@ -95,8 +104,9 @@ class WentzelDistribution // Sample the scattered polar angle template - inline CELER_FUNCTION real_type sample_cos_t(real_type cos_t_max, - Engine& rng) const; + inline CELER_FUNCTION real_type sample_costheta(real_type cos_thetamin, + real_type cos_thetamax, + Engine& rng) const; // Helper function for calculating the flat form factor inline static CELER_FUNCTION real_type flat_form_factor(real_type x); @@ -113,17 +123,25 @@ class WentzelDistribution */ CELER_FUNCTION WentzelDistribution::WentzelDistribution( + NativeCRef const& wentzel, + WentzelHelper const& helper, ParticleTrackView const& particle, IsotopeView const& target, - CoulombScatteringElementData const& element_data, - Energy cutoff, - CoulombScatteringRef const& data) - : data_(data) + ElementId el_id, + real_type cos_thetamin, + real_type cos_thetamax) + : wentzel_(wentzel) + , helper_(helper) , particle_(particle) , target_(target) - , element_data_(element_data) - , helper_(particle, target.atomic_number(), data, cutoff) + , el_id_(el_id) + , cos_thetamin_(cos_thetamin) + , cos_thetamax_(cos_thetamax) { + CELER_EXPECT(el_id_ < wentzel_.elem_data.size()); + CELER_EXPECT(cos_thetamin_ >= -1 && cos_thetamin_ <= 1); + CELER_EXPECT(cos_thetamax_ >= -1 && cos_thetamax_ <= 1); + CELER_EXPECT(cos_thetamax_ <= cos_thetamin_); } //---------------------------------------------------------------------------// @@ -134,25 +152,28 @@ template CELER_FUNCTION real_type WentzelDistribution::operator()(Engine& rng) const { real_type cos_theta = 1; - if (BernoulliDistribution(helper_.calc_xs_ratio())(rng)) + if (BernoulliDistribution( + helper_.calc_xs_ratio(cos_thetamin_, cos_thetamax_))(rng)) { // Scattered off of electrons - cos_theta = this->sample_cos_t(helper_.costheta_max_electron(), rng); + real_type const cos_thetamax_elec = helper_.cos_thetamax_electron(); + real_type cos_thetamin = max(cos_thetamin_, cos_thetamax_elec); + real_type cos_thetamax = max(cos_thetamax_, cos_thetamax_elec); + CELER_ASSERT(cos_thetamin > cos_thetamax); + cos_theta = this->sample_costheta(cos_thetamin, cos_thetamax, rng); } else { // Scattered off of nucleus - cos_theta = this->sample_cos_t(-1, rng); + cos_theta = this->sample_costheta(cos_thetamin_, cos_thetamax_, rng); // Calculate rejection for fake scattering // TODO: Reference? - real_type mott_coeff - = 1 + real_type(2e-4) * ipow<2>(target_.atomic_number().get()); - MottRatioCalculator mott_xsec(element_data_, + MottRatioCalculator mott_xsec(wentzel_.elem_data[el_id_], std::sqrt(particle_.beta_sq())); real_type g_rej = mott_xsec(cos_theta) * ipow<2>(this->calculate_form_factor(cos_theta)) - / mott_coeff; + / helper_.mott_factor(); if (!BernoulliDistribution(g_rej)(rng)) { @@ -160,7 +181,6 @@ CELER_FUNCTION real_type WentzelDistribution::operator()(Engine& rng) const cos_theta = 1; } } - return cos_theta; } @@ -179,7 +199,7 @@ CELER_FUNCTION real_type WentzelDistribution::calculate_form_factor(real_type cos_t) const { real_type mt_sq = this->mom_transfer_sq(cos_t); - switch (data_.form_factor_type) + switch (wentzel_.params.form_factor_type) { case NuclearFormFactorType::none: return 1; @@ -194,8 +214,9 @@ WentzelDistribution::calculate_form_factor(real_type cos_t) const return 1 / ipow<2>(1 + this->nuclear_form_prefactor() * mt_sq); case NuclearFormFactorType::gaussian: return std::exp(-2 * this->nuclear_form_prefactor() * mt_sq); + default: + CELER_ASSERT_UNREACHABLE(); } - CELER_ASSERT_UNREACHABLE(); } //---------------------------------------------------------------------------// @@ -241,28 +262,36 @@ CELER_FUNCTION real_type WentzelDistribution::mom_transfer_sq(real_type cos_t) c */ CELER_FUNCTION real_type WentzelDistribution::nuclear_form_prefactor() const { - return data_.nuclear_form_prefactor[target_.isotope_id()]; + CELER_EXPECT(target_.isotope_id() < wentzel_.nuclear_form_prefactor.size()); + return wentzel_.nuclear_form_prefactor[target_.isotope_id()]; } //---------------------------------------------------------------------------// /*! * Sample the scattering polar angle of the incident particle. * - * The probability is given in [Fern] - * eqn 88 and is nomalized on the interval - * \f$ cos\theta \in [1, \cos\theta_{max}] \f$. + * The probability is given in [Fern] eqn 88 and is nomalized on the interval + * \f$ cos\theta \in [\cos\theta_{min}, \cos\theta_{max}] \f$. The sampling + * function for \f$ \mu = \frac{1}{2}(1 - \cos\theta) \f$ is + * \f[ + \mu = \mu_1 + \frac{(A + \mu_1) \xi (\mu_2 - \mu_1)}{A + \mu_2 - \xi (\mu_2 + - \mu_1)}, + * \f] + * where \f$ \mu_1 = \frac{1}{2}(1 - \cos\theta_{min}) \f$, \f$ \mu_2 = + * \frac{1}{2}(1 - \cos\theta_{max}) \f$, \f$ A \f$ is the screening + * coefficient, and \f$ \xi \sim U(0,1) \f$. */ template -CELER_FUNCTION real_type WentzelDistribution::sample_cos_t(real_type cos_t_max, - Engine& rng) const +CELER_FUNCTION real_type WentzelDistribution::sample_costheta( + real_type cos_thetamin, real_type cos_thetamax, Engine& rng) const { // Sample scattering angle [Fern] eqn 92, where cos(theta) = 1 - 2*mu - // For incident electrons / positrons, theta_min = 0 always - real_type const mu = real_type{0.5} * (1 - cos_t_max); - real_type const xi = generate_canonical(rng); + real_type const mu1 = real_type{0.5} * (1 - cos_thetamin); + real_type const mu2 = real_type{0.5} * (1 - cos_thetamax); + real_type const w = generate_canonical(rng) * (mu2 - mu1); real_type const sc = helper_.screening_coefficient(); - real_type result = 1 - 2 * sc * mu * (1 - xi) / (sc + mu * xi); + real_type result = 1 - 2 * mu1 - 2 * (sc + mu1) * w / (sc + mu2 - w); CELER_ENSURE(result >= -1 && result <= 1); return result; } diff --git a/src/celeritas/em/executor/CoulombScatteringExecutor.hh b/src/celeritas/em/executor/CoulombScatteringExecutor.hh index 85895ea6fc..81b3f9b6e0 100644 --- a/src/celeritas/em/executor/CoulombScatteringExecutor.hh +++ b/src/celeritas/em/executor/CoulombScatteringExecutor.hh @@ -10,6 +10,7 @@ #include "corecel/Assert.hh" #include "corecel/Macros.hh" #include "celeritas/em/data/CoulombScatteringData.hh" +#include "celeritas/em/data/WentzelOKVIData.hh" #include "celeritas/em/interactor/CoulombScatteringInteractor.hh" #include "celeritas/global/CoreTrackView.hh" #include "celeritas/mat/IsotopeSelector.hh" @@ -27,7 +28,8 @@ struct CoulombScatteringExecutor inline CELER_FUNCTION Interaction operator()(celeritas::CoreTrackView const& track); - CoulombScatteringRef params; + CoulombScatteringData params; + NativeCRef wentzel; }; //---------------------------------------------------------------------------// @@ -56,7 +58,7 @@ CoulombScatteringExecutor::operator()(CoreTrackView const& track) // Construct the interactor CoulombScatteringInteractor interact( - params, particle, dir, target, element_id, cutoffs); + params, wentzel, particle, dir, material, target, element_id, cutoffs); // Execute the interactor return interact(rng); diff --git a/src/celeritas/em/interactor/CoulombScatteringInteractor.hh b/src/celeritas/em/interactor/CoulombScatteringInteractor.hh index 6b7cc9d582..89929a9ea9 100644 --- a/src/celeritas/em/interactor/CoulombScatteringInteractor.hh +++ b/src/celeritas/em/interactor/CoulombScatteringInteractor.hh @@ -13,6 +13,7 @@ #include "celeritas/Quantities.hh" #include "celeritas/Types.hh" #include "celeritas/em/data/CoulombScatteringData.hh" +#include "celeritas/em/data/WentzelOKVIData.hh" #include "celeritas/em/distribution/WentzelDistribution.hh" #include "celeritas/mat/ElementView.hh" #include "celeritas/mat/MaterialView.hh" @@ -52,11 +53,13 @@ class CoulombScatteringInteractor public: //! Construct with shared and state data inline CELER_FUNCTION - CoulombScatteringInteractor(CoulombScatteringRef const& shared, + CoulombScatteringInteractor(CoulombScatteringData const& shared, + NativeCRef const& wentzel, ParticleTrackView const& particle, Real3 const& inc_direction, + MaterialView const& material, IsotopeView const& target, - ElementId const& el_id, + ElementId el_id, CutoffView const& cutoffs); //! Sample an interaction with the given RNG @@ -75,8 +78,11 @@ class CoulombScatteringInteractor // Target isotope IsotopeView const& target_; + // Helper for calculating xs ratio and other quantities + WentzelHelper const helper_; + // Scattering direction distribution of the Wentzel model - WentzelDistribution const sample_angle; + WentzelDistribution const sample_angle_; //// HELPER FUNCTIONS //// @@ -92,21 +98,32 @@ class CoulombScatteringInteractor */ CELER_FUNCTION CoulombScatteringInteractor::CoulombScatteringInteractor( - CoulombScatteringRef const& shared, + CoulombScatteringData const& shared, + NativeCRef const& wentzel, ParticleTrackView const& particle, Real3 const& inc_direction, + MaterialView const& material, IsotopeView const& target, - ElementId const& el_id, + ElementId el_id, CutoffView const& cutoffs) : inc_direction_(inc_direction) , particle_(particle) , target_(target) - , sample_angle(particle, - target, - shared.elem_data[el_id], - // TODO: Use proton when supported - cutoffs.energy(shared.ids.electron), - shared) + , helper_(particle, + material, + target.atomic_number(), + wentzel, + shared.ids, + // TODO: Use the proton production cutoff when the recoiled + // nucleus production is supported + cutoffs.energy(shared.ids.electron)) + , sample_angle_(wentzel, + helper_, + particle, + target, + el_id, + helper_.cos_thetamax_nuclear(), + shared.cos_thetamax()) { CELER_EXPECT(particle_.particle_id() == shared.ids.electron || particle_.particle_id() == shared.ids.positron); @@ -125,7 +142,7 @@ CELER_FUNCTION Interaction CoulombScatteringInteractor::operator()(Engine& rng) Interaction result; // Sample the new direction - real_type cos_theta = sample_angle(rng); + real_type cos_theta = sample_angle_(rng); UniformRealDistribution sample_phi(0, 2 * constants::pi); result.direction = rotate(inc_direction_, from_spherical(cos_theta, sample_phi(rng))); diff --git a/src/celeritas/em/model/CoulombScatteringModel.cc b/src/celeritas/em/model/CoulombScatteringModel.cc index 59d2219072..9da149e6dc 100644 --- a/src/celeritas/em/model/CoulombScatteringModel.cc +++ b/src/celeritas/em/model/CoulombScatteringModel.cc @@ -7,35 +7,32 @@ //---------------------------------------------------------------------------// #include "CoulombScatteringModel.hh" +#include + #include "celeritas_config.h" -#include "corecel/sys/ScopedMem.hh" +#include "celeritas/Constants.hh" #include "celeritas/Units.hh" #include "celeritas/em/data/CoulombScatteringData.hh" #include "celeritas/em/executor/CoulombScatteringExecutor.hh" #include "celeritas/em/interactor/detail/PhysicsConstants.hh" -#include "celeritas/em/process/CoulombScatteringProcess.hh" +#include "celeritas/em/params/WentzelOKVIParams.hh" #include "celeritas/global/ActionLauncher.hh" #include "celeritas/global/CoreParams.hh" #include "celeritas/global/TrackExecutor.hh" #include "celeritas/io/ImportParameters.hh" #include "celeritas/io/ImportProcess.hh" -#include "celeritas/mat/ElementView.hh" -#include "celeritas/mat/MaterialParams.hh" #include "celeritas/phys/InteractionApplier.hh" #include "celeritas/phys/PDGNumber.hh" #include "celeritas/phys/ParticleParams.hh" -#include "celeritas/phys/ParticleView.hh" namespace celeritas { //---------------------------------------------------------------------------// /*! - * Construct from model ID, particle data, and material data. + * Construct from model ID and shared data. */ CoulombScatteringModel::CoulombScatteringModel(ActionId id, ParticleParams const& particles, - MaterialParams const& materials, - Options const& options, SPConstImported data) : imported_(data, particles, @@ -45,32 +42,16 @@ CoulombScatteringModel::CoulombScatteringModel(ActionId id, { CELER_EXPECT(id); - ScopedMem record_mem("CoulombScatteringModel.construct"); - HostVal host_data; - - // This is where the data is built and transfered to the device - host_data.ids.action = id; - host_data.ids.electron = particles.find(pdg::electron()); - host_data.ids.positron = particles.find(pdg::positron()); - host_data.ids.proton = particles.find(pdg::proton()); - - CELER_VALIDATE(host_data.ids, - << "missing IDs (required for " << this->description() - << ")"); - - // Select form factor - host_data.form_factor_type = options.form_factor_model; + data_.action = id; + data_.ids.electron = particles.find(pdg::electron()); + data_.ids.positron = particles.find(pdg::positron()); - // Pass user-defined screening factor - host_data.screening_factor = options.screening_factor; + CELER_VALIDATE(data_.ids, + << "missing electron and/or positron particles (required " + "for " + << this->description() << ")"); - // Load Mott coefficients - build_data(host_data, materials); - - // Construct data on device - data_ = CollectionMirror{std::move(host_data)}; - - CELER_ENSURE(this->data_); + CELER_ENSURE(data_); } //---------------------------------------------------------------------------// @@ -110,11 +91,14 @@ auto CoulombScatteringModel::micro_xs(Applicability applic) const void CoulombScatteringModel::execute(CoreParams const& params, CoreStateHost& state) const { + CELER_EXPECT(params.wentzel()); + auto execute = make_action_track_executor( params.ptr(), state.ptr(), this->action_id(), - InteractionApplier{CoulombScatteringExecutor{this->host_ref()}}); + InteractionApplier{CoulombScatteringExecutor{ + this->host_ref(), params.wentzel()->host_ref()}}); return launch_action(*this, params, state, execute); } @@ -133,842 +117,7 @@ void CoulombScatteringModel::execute(CoreParams const&, CoreStateDevice&) const */ ActionId CoulombScatteringModel::action_id() const { - return this->host_ref().ids.action; -} - -//---------------------------------------------------------------------------// -/*! - * Load Mott coefficients and construct per-element data. - */ -void CoulombScatteringModel::build_data( - HostVal& host_data, MaterialParams const& materials) -{ - // Build element data - size_type const num_elements = materials.num_elements(); - auto elem_data = make_builder(&host_data.elem_data); - elem_data.reserve(num_elements); - - for (auto el_id : range(ElementId{num_elements})) - { - // Load Mott coefficients - CoulombScatteringElementData z_data; - z_data.mott_coeff - = get_mott_coeff_matrix(materials.get(el_id).atomic_number()); - elem_data.push_back(z_data); - } - - auto prefactors = make_builder(&host_data.nuclear_form_prefactor); - prefactors.reserve(materials.num_isotopes()); - for (auto iso_id : range(IsotopeId{materials.num_isotopes()})) - { - prefactors.push_back( - this->calc_nuclear_form_prefactor(materials.get(iso_id))); - } -} - -//---------------------------------------------------------------------------// -/*! - * Get interpolated Mott coefficients for an element. - * - * These coefficients are used by the Lijian, Quing, Zhengming - * expression [PRM 8.48] for atomic numbers 1 <= Z <= 92. - * This data was taken from Geant4's G4MottData.hh file. - * - * For higher Z values, the PRM cites a numerical solution by Boschini - * et al (2013), but does not implement it. - */ -CoulombScatteringElementData::MottCoeffMatrix -CoulombScatteringModel::get_mott_coeff_matrix(AtomicNumber z) -{ - CELER_EXPECT(z); - - using CoeffMat = CoulombScatteringElementData::MottCoeffMatrix; - // clang-format off - static CoeffMat const mott_coeffs[] = { - //H....................................... - CoeffMat{{ - {1.e0,2.67363e-8,7.1153e-8,-9.7703e-8,-6.69132e-7,-3.09263e-7}, - {1.17182e-2,1.62222e-2,-5.90397e-5,-1.05585e-4,4.17873e-4,9.13843e-4}, - {-2.65955e-1,-7.29531e-1,-4.99796e-1,2.83507e-4,-9.09042e-4,-2.20244e-3}, - {-1.82348e-4,-8.86355e-5,-1.90554e-4,-2.49708e-4,6.35004e-4,1.73523e-3}, - {4.70966e-5,-4.09705e-6,3.75218e-5,8.05645e-5,-1.90534e-4,-5.42847e-4} - }}, - //He....................................... - CoeffMat{{ - {1.e0,3.76476e-8,-3.05313e-7,-3.27422e-7,2.44235e-6,4.08754e-6}, - {2.35767e-2,3.24642e-2,-6.37269e-4,-7.6916e-4,5.28004e-3,9.45642e-3}, - {-2.73743e-1,-7.40767e-1,-4.98195e-1,1.74337e-3,-1.25798e-2,-2.24046e-2}, - {-7.79128e-4,-4.14495e-4,-1.62657e-3,-1.37286e-3,1.04319e-2,1.83488e-2}, - {2.02855e-4,1.94598e-6,4.30102e-4,4.3218e-4,-3.31526e-3,-5.81788e-3} - }}, - //Li.................................................. - CoeffMat{{ - {1.e0,7.00357e-8,-3.15076e-7,-4.24915e-7,2.45516e-6,4.90187e-6}, - {3.55657e-2,4.87956e-2,-1.95525e-3,-2.7866e-3,1.6549e-2,3.11496e-2}, - {-2.81171e-1,-7.52015e-1,-4.95329e-1,5.83548e-3,-3.3983e-2,-6.55379e-2}, - {-1.83452e-3,-8.12746e-4,-3.84675e-3,-4.44467e-3,2.55871e-2,4.99483e-2}, - {4.79031e-4,-3.89615e-5,1.01022e-3,1.39133e-3,-7.99398e-3,-1.56366e-2} - }}, - //Be.................................................. - CoeffMat{{ - {1.e0,7.58881e-8,4.705e-8,2.48041e-7,-2.06053e-6,-1.97319e-6}, - {4.76788e-2,6.522e-2,-4.54331e-3,-6.50318e-3,3.76564e-2,7.17176e-2}, - {-2.88203e-1,-7.63217e-1,-4.90337e-1,1.22839e-2,-6.86398e-2,-1.35769e-1}, - {-3.37733e-3,-1.36514e-3,-7.51614e-3,-8.78592e-3,4.78572e-2,9.69021e-2}, - {8.81822e-4,-1.02577e-4,1.99797e-3,2.72661e-3,-1.48296e-2,-3.0106e-2} - }}, - //B.................................................... - CoeffMat{{ - {9.99999e-1,7.91498e-8,1.84164e-6,2.68534e-6,-1.8163e-5,-2.69021e-5}, - {5.98818e-2,8.17654e-2,-7.70811e-3,-1.12378e-2,6.38329e-2,1.25339e-1}, - {-2.94716e-1,-7.74405e-1,-4.8622e-1,1.77367e-2,-9.46825e-2,-2.01789e-1}, - {-5.52375e-3,-2.05348e-3,-9.44915e-3,-1.08135e-2,5.41024e-2,1.25257e-1}, - {1.44555e-3,-1.99404e-4,2.36742e-3,3.29655e-3,-1.64122e-2,-3.8375e-2} - }}, - //C................................................... - CoeffMat{{ - {9.99999e-1,7.68158e-8,5.18185e-6,7.34245e-6,-4.9478e-5,-7.71923e-5}, - {7.21461e-2,9.84618e-2,-1.06535e-2,-1.62358e-2,8.59238e-2,1.78727e-1}, - {-3.00622e-1,-7.85616e-1,-4.85735e-1,1.91563e-2,-8.10204e-2,-2.15386e-1}, - {-8.34809e-3,-2.85241e-3,-7.03252e-3,-7.56786e-3,1.44975e-2,8.79093e-2}, - {2.18964e-3,-3.42022e-4,1.3293e-3,2.20108e-3,-3.57927e-3,-2.5928e-2} - }}, - //N.................................................... - CoeffMat{{ - {9.99999e-1,8.36312e-8,1.09116e-5,1.47812e-5,-1.02733e-4,-1.62724e-4}, - {8.44142e-2,1.1531e-1,-1.1723e-2,-1.94732e-2,8.92604e-2,2.09303e-1}, - {-3.05743e-1,-7.96809e-1,-4.93957e-1,1.01607e-2,1.67761e-2,-1.05909e-1}, - {-1.2009e-2,-3.80678e-3,4.51195e-3,6.93472e-3,-1.12405e-1,-8.15484e-2}, - {3.16048e-3,-5.22237e-4,-2.58261e-3,-2.38303e-3,3.63393e-2,2.75127e-2} - }}, - //O................................................... - CoeffMat{{ - {9.99998e-1,1.57323e-8,1.77595e-5,2.56082e-5,-1.67537e-4,-2.73755e-4}, - {9.66438e-2,1.32264e-1,-9.53841e-3,-1.83707e-2,6.01664e-2,1.93357e-1}, - {-3.09969e-1,-8.0779e-1,-5.14392e-1,-1.67153e-2,2.3387e-1,1.916e-1}, - {-1.65906e-2,-5.11585e-3,2.80424e-2,3.94663e-2,-3.5572e-1,-4.39251e-1}, - {4.37866e-3,-6.81795e-4,-1.0152e-2,-1.24875e-2,1.11484e-1,1.38105e-1} - }}, - //F................................................. - CoeffMat{{ - {9.99997e-1,-8.06132e-8,2.49797e-5,3.8512e-5,-2.37451e-4,-3.99607e-4}, - {1.08782e-1,1.49306e-1,-2.50975e-3,-1.05471e-2,-1.64831e-2,1.05733e-1}, - {-3.13165e-1,-8.18489e-1,-5.50832e-1,-6.74447e-2,6.06357e-1,7.39717e-1}, - {-2.21976e-2,-6.84023e-3,6.65411e-2,9.48702e-2,-7.43989e-1,-1.03582e0}, - {5.87088e-3,-8.1048e-4,-2.21731e-2,-2.94422e-2,2.2954e-1,3.19669e-1} - }}, - //Ne................................................. - CoeffMat{{ - {9.99997e-1,-1.87404e-7,3.10276e-5,5.20007e-5,-2.98132e-4,-5.19259e-4}, - {1.20783e-1,1.66407e-1,1.06608e-2,6.48772e-3,-1.53031e-1,-7.59354e-2}, - {-3.15222e-1,-8.28793e-1,-6.0574e-1,-1.47812e-1,1.1576e0,1.58565e0}, - {-2.89055e-2,-9.08096e-3,1.21467e-1,1.77575e-1,-1.2911e0,-1.90333e0}, - {7.65342e-3,-8.85417e-4,-3.89092e-2,-5.4404e-2,3.93087e-1,5.79439e-1} - }}, - //Na................................................. - CoeffMat{{ - {9.99996e-1,-2.44548e-7,3.31019e-5,6.29483e-5,-3.24667e-4,-5.95527e-4}, - {1.32615e-1,1.83566e-1,3.04158e-2,3.40925e-2,-3.54681e-1,-3.63044e-1}, - {-3.16092e-1,-8.38704e-1,-6.78558e-1,-2.59346e-1,1.88547e0,2.73632e0}, - {-3.67233e-2,-1.18139e-2,1.91089e-1,2.87408e-1,-1.98397e0,-3.03075e0}, - {9.72033e-3,-9.2638e-4,-5.95654e-2,-8.69829e-2,5.95744e-1,9.10242e-1} - }}, - //Mg................................................. - CoeffMat{{ - {9.99995e-1,-2.12227e-7,2.95645e-5,6.92848e-5,-3.02153e-4,-6.05145e-4}, - {1.44258e-1,2.00775e-1,5.67845e-2,7.35166e-2,-6.22861e-1,-7.62213e-1}, - {-3.15754e-1,-8.48196e-1,-7.67318e-1,-4.02984e-1,2.77477e0,4.18114e0}, - {-4.56307e-2,-1.50425e-2,2.72232e-1,4.23528e-1,-2.79606e0,-4.38863e0}, - {1.2056e-2,-9.44637e-4,-8.28738e-2,-1.26564e-1,8.26726e-1,1.29882e0} - }}, - //Al..................................................... - CoeffMat{{ - {9.99995e-1,-4.03407e-8,1.86047e-5,6.85201e-5,-2.14503e-4,-5.22528e-4}, - {1.55704e-1,2.18048e-1,8.88994e-2,1.24878e-1,-9.51331e-1,-1.26824e0}, - {-3.14244e-1,-8.57322e-1,-8.6719e-1,-5.75787e-1,3.78571e0,5.87052e0}, - {-5.55526e-2,-1.86861e-2,3.5886e-1,5.81094e-1,-3.67623e0,-5.908e0}, - {1.46269e-2,-9.79742e-4,-1.06652e-1,-1.71226e-1,1.06737e0,1.71918e0} - }}, - //Si..................................................... - CoeffMat{{ - {9.99994e-1,3.00267e-7,-1.1184e-6,5.88256e-5,-4.78456e-5,-3.25731e-4}, - {1.6696e-1,2.35405e-1,1.25215e-1,1.87646e-1,-1.32685e0,-1.86549e0}, - {-3.1163e-1,-8.66152e-1,-9.71254e-1,-7.72715e-1,4.85654e0,7.7215e0}, - {-6.63778e-2,-2.26481e-2,4.42898e-1,7.53182e-1,-4.55172e0,-7.4867e0}, - {1.73883e-2,-1.07669e-3,-1.28075e-1,-2.18389e-1,1.29217e0,2.13475e0} - }}, - //P..................................................... - CoeffMat{{ - {9.99994e-1,8.94829e-7,-2.98434e-5,3.82193e-5,2.00584e-4,-6.40482e-6}, - {1.78039e-1,2.52912e-1,1.63761e-1,2.60132e-1,-1.73287e0,-2.53185e0}, - {-3.08007e-1,-8.74905e-1,-1.07146e0,-9.85062e-1,5.91697e0,9.63265e0}, - {-7.79747e-2,-2.66797e-2,5.15288e-1,9.29261e-1,-5.34252e0,-9.00574e0}, - {2.02892e-2,-1.33011e-3,-1.44039e-1,-2.6433e-1,1.4736e0,2.50398e0} - }}, - //S.................................................... - CoeffMat{{ - {9.99994e-1,1.75397e-6,-6.7331e-5,6.29524e-6,5.29623e-4,4.35288e-4}, - {1.88968e-1,2.70612e-1,2.01975e-1,3.40574e-1,-2.14737e0,-3.23836e0}, - {-3.03499e-1,-8.83717e-1,-1.15816e0,-1.20414e0,6.88176e0,1.14841e1}, - {-9.01806e-2,-3.06202e-2,5.65581e-1,1.09902e0,-5.95552e0,-1.03302e1}, - {2.32694e-2,-1.80614e-3,-1.51041e-1,-3.05449e-1,1.58037e0,2.78083e0} - }}, - //Cl.................................................... - CoeffMat{{ - {9.99994e-1,3.07931e-6,-1.11876e-4,-4.10164e-5,9.17095e-4,9.80145e-4}, - {1.99765e-1,2.88611e-1,2.37501e-1,4.25803e-1,-2.55105e0,-3.95585e0}, - {-2.98206e-1,-8.9293e-1,-1.22279e0,-1.4169e0,7.67836e0,1.31601e1}, - {-1.02865e-1,-3.40967e-2,5.84677e-1,1.24786e0,-6.31301e0,-1.13328e1}, - {2.628e-2,-2.63995e-3,-1.46076e-1,-3.36795e-1,1.58677e0,2.92251e0} - }}, - //Ar.................................................... - CoeffMat{{ - {9.99993e-1,4.49776e-6,-1.65136e-4,-9.76754e-5,1.39664e-3,1.66293e-3}, - {2.10469e-1,3.06924e-1,2.66793e-1,5.13797e-1,-2.90958e0,-4.63816e0}, - {-2.92294e-1,-9.0256e-1,-1.25307e0,-1.6147e0,8.18574e0,1.44912e1}, - {-1.15831e-1,-3.70891e-2,5.59807e-1,1.36619e0,-6.28824e0,-1.18327e1}, - {2.92513e-2,-3.84903e-3,-1.24976e-1,-3.55149e-1,1.45127e0,2.86925e0} - }}, - //K................................................. - CoeffMat{{ - {9.99993e-1,6.01488e-6,-2.22125e-4,-1.61774e-4,1.92058e-3,2.41975e-3}, - {2.21091e-1,3.25648e-1,2.87732e-1,6.01632e-1,-3.20436e0,-5.25724e0}, - {-2.85814e-1,-9.12889e-1,-1.24213e0,-1.78635e0,8.34196e0,1.53776e1}, - {-1.29005e-1,-3.92986e-2,4.84255e-1,1.44206e0,-5.81999e0,-1.17275e1}, - {3.21555e-2,-5.54272e-3,-8.56572e-2,-3.56589e-1,1.1548e0,2.58829e0} - }}, - //Ca................................................. - CoeffMat{{ - {9.99993e-1,8.01467e-6,-2.79242e-4,-2.3682e-4,2.45459e-3,3.21683e-3}, - {2.31651e-1,3.44948e-1,2.9782e-1,6.85187e-1,-3.41294e0,-5.77715e0}, - {-2.78858e-1,-9.24428e-1,-1.18215e0,-1.91691e0,8.07489e0,1.56969e1}, - {-1.42276e-1,-4.01888e-2,3.50466e-1,1.45983e0,-4.83806e0,-1.08936e1}, - {3.49529e-2,-7.90933e-3,-2.58002e-2,-3.36028e-1,6.7574e-1,2.04052e0} - }}, - //Sc................................................. - CoeffMat{{ - {9.99992e-1,1.04277e-5,-3.35126e-4,-3.21042e-4,2.98507e-3,4.03325e-3}, - {2.42172e-1,3.64954e-1,2.94606e-1,7.60693e-1,-3.51409e0,-6.1646e0}, - {-2.71512e-1,-9.37543e-1,-1.0657e0,-1.99328e0,7.31863e0,1.53396e1}, - {-1.5554e-1,-3.93862e-2,1.51477e-1,1.40614e0,-3.28024e0,-9.22338e0}, - {3.76066e-2,-1.10812e-2,5.66831e-2,-2.8921e-1,-4.60274e-3,1.19273e0} - }}, - //Ti................................................. - CoeffMat{{ - {9.99992e-1,1.30838e-5,-3.8407e-4,-4.09294e-4,3.47025e-3,4.81071e-3}, - {2.52646e-1,3.85718e-1,2.7633e-1,8.25665e-1,-3.48888e0,-6.39017e0}, - {-2.63758e-1,-9.52326e-1,-8.88179e-1,-2.0072e0,6.0196e0,1.42157e1}, - {-1.68806e-1,-3.68095e-2,-1.16429e-1,1.27309e0,-1.09991e0,-6.63409e0}, - {4.01184e-2,-1.50938e-2,1.6274e-1,-2.1376e-1,-8.99142e-1,2.08129e-2} - }}, - //V................................................. - CoeffMat{{ - {9.99991e-1,1.59363e-5,-4.27315e-4,-5.01461e-4,3.91365e-3,5.55096e-3}, - {2.63096e-1,4.07357e-1,2.40649e-1,8.7645e-1,-3.31641e0,-6.42069e0}, - {-2.55682e-1,-9.6909e-1,-6.43149e-1,-1.94678e0,4.11915e0,1.22255e1}, - {-1.81974e-1,-3.21462e-2,-4.58817e-1,1.04913e0,1.75388e0,-3.03434e0}, - {4.24541e-2,-2.00595e-2,2.93916e-1,-1.06138e-1,-2.02203e0,-1.50181e0} - }}, - //Cr................................................. - CoeffMat{{ - {9.9999e-1,1.8895e-5,-4.59994e-4,-5.93663e-4,4.27684e-3,6.2009e-3}, - {2.73504e-1,4.2999e-1,1.86473e-1,9.09921e-1,-2.98441e0,-6.23344e0}, - {-2.47225e-1,-9.88118e-1,-3.28759e-1,-1.80252e0,1.5909e0,9.30968e0}, - {-1.951e-1,-2.51242e-2,-8.76244e-1,7.25592e-1,5.2965e0,1.62177e0}, - {4.46307e-2,-2.60754e-2,4.50054e-1,3.61638e-2,-3.37524e0,-3.38601e0} - }}, - //Mn................................................. - CoeffMat{{ - {9.99989e-1,2.18906e-5,-4.79199e-4,-6.83164e-4,4.53467e-3,6.72485e-3}, - {2.83854e-1,4.53722e-1,1.12877e-1,9.23203e-1,-2.48247e0,-5.80855e0}, - {-2.38338e-1,-1.00965e0,5.61453e-2,-1.56605e0,-1.58353e0,5.42138e0}, - {-2.08226e-1,-1.55138e-2,-1.36846e0,2.95168e-1,9.53392e0,7.3655e0}, - {4.66619e-2,-3.32255e-2,6.30711e-1,2.15167e-1,-4.95748e0,-5.63744e0} - }}, - //Fe................................................... - CoeffMat{{ - {9.99987e-1,2.482e-5,-4.82895e-4,-7.67488e-4,4.669e-3,7.09581e-3}, - {2.94123e-1,4.78653e-1,1.93256e-2,9.13569e-1,-1.80354e0,-5.1304e0}, - {-2.28945e-1,-1.0339e0,5.1121e-1,-1.22996e0,-5.40939e0,5.3149e-1}, - {-2.21426e-1,-3.12679e-3,-1.93353e0,-2.48179e-1,1.44576e1,1.42077e1}, - {4.85718e-2,-4.15796e-2,8.34879e-1,4.32418e-1,-6.76244e0,-8.25464e0} - }}, - //Co.................................................. - CoeffMat{{ - {9.99985e-1,2.76168e-5,-4.67522e-4,-8.43978e-4,4.65008e-3,7.27476e-3}, - {3.0428e-1,5.04856e-1,-9.42256e-2,8.78905e-1,-9.44256e-1,-4.18847e0}, - {-2.18939e-1,-1.06098e0,1.03431e0,-7.89123e-1,-9.87815e0,-5.36997e0}, - {-2.34805e-1,1.21281e-2,-2.56765e0,-9.08015e-1,2.00436e1,2.21379e1}, - {5.03942e-2,-5.11768e-2,1.0609e0,6.88631e-1,-8.77867e0,-1.12289e1} - }}, - //Ni.................................................. - CoeffMat{{ - {9.99982e-1,3.00792e-5,-4.33447e-4,-9.09366e-4,4.47786e-3,7.25072e-3}, - {3.14283e-1,5.32444e-1,-2.27528e-1,8.16951e-1,9.53704e-2,-2.9773e0}, - {-2.08186e-1,-1.09112e0,1.62237e0,-2.38394e-1,-1.49699e1,-1.22743e1}, - {-2.48493e-1,3.04543e-2,-3.26601e0,-1.6876e0,2.62568e1,3.11259e1}, - {5.21709e-2,-6.20908e-2,1.30683e0,9.84353e-1,-1.09909e1,-1.45449e1} - }}, - //Cu................................................... - CoeffMat{{ - {9.99979e-1,3.24569e-5,-3.76717e-4,-9.64e-4,4.11997e-3,6.98761e-3}, - {3.24082e-1,5.61534e-1,-3.79993e-1,7.25313e-1,1.31342e0,-1.49163e0}, - {-1.96521e-1,-1.12457e0,2.27091e0,4.27678e-1,-2.06558e1,-2.0169e1}, - {-2.62655e-1,5.20812e-2,-4.02224e0,-2.59048e0,3.30508e1,4.1135e1}, - {5.39556e-2,-7.44077e-2,1.57015e0,1.32022e0,-1.33802e1,-1.81848e1} - }}, - //Zn................................................... - CoeffMat{{ - {9.99976e-1,3.39628e-5,-2.98845e-4,-9.99651e-4,3.58523e-3,6.47782e-3}, - {3.33647e-1,5.92068e-1,-5.5102e-1,6.0363e-1,2.70712e0,2.67853e-1}, - {-1.83849e-1,-1.16095e0,2.97559e0,1.20728e0,-2.69053e1,-2.90214e1}, - {-2.77375e-1,7.65838e-2,-4.83023e0,-3.61238e0,4.03786e1,5.21084e1}, - {5.57745e-2,-8.79952e-2,1.84848e0,1.69425e0,-1.59274e1,-2.21242e1} - }}, - //Ga................................................... - CoeffMat{{ - {9.99972e-1,3.48473e-5,-1.97828e-4,-1.01659e-3,2.85509e-3,5.69725e-3}, - {3.42904e-1,6.24164e-1,-7.39023e-1,4.50306e-1,4.26553e0,2.29299e0}, - {-1.6994e-1,-1.20051e0,3.72868e0,2.10276e0,-3.36594e1,-3.87726e1}, - {-2.92882e-1,1.04194e-1,-5.68033e0,-4.75339e0,4.81633e1,6.39614e1}, - {5.77008e-2,-1.02941e-1,2.13826e0,2.10589e0,-1.86034e1,-2.63294e1} - }}, - //Ge................................................. - CoeffMat{{ - {9.99968e-1,3.47804e-5,-7.11898e-5,-1.01028e-3,1.90919e-3,4.61426e-3}, - {3.51793e-1,6.57815e-1,-9.42175e-1,2.6494e-1,5.97606e0,4.57222e0}, - {-1.54595e-1,-1.24305e0,4.5217e0,3.11205e0,-4.08539e1,-4.93518e1}, - {-3.09367e-1,1.34661e-1,-6.56214e0,-6.00882e0,5.63228e1,7.65965e1}, - {5.97948e-2,-1.19174e-1,2.4357e0,2.553e0,-2.13779e1,-3.07628e1} - }}, - //As................................................. - CoeffMat{{ - {9.99963e-1,3.37519e-5,8.13037e-5,-9.78638e-4,7.41412e-4,3.21498e-3}, - {3.60246e-1,6.93065e-1,-1.1588e0,4.68519e-2,7.82662e0,7.09456e0}, - {-1.37615e-1,-1.28855e0,5.34673e0,4.23402e0,-4.84263e1,-6.06897e1}, - {-3.27017e-1,1.67941e-1,-7.46593e0,-7.37492e0,6.47773e1,8.99181e1}, - {6.21158e-2,-1.36692e-1,2.73726e0,3.03374e0,-2.42209e1,-3.53873e1} - }}, - //Se................................................. - CoeffMat{{ - {9.99958e-1,3.14983e-5,2.59741e-4,-9.19008e-4,-6.49202e-4,1.49008e-3}, - {3.68196e-1,7.29888e-1,-1.38664e0,-2.03807e-1,9.80062e0,9.84155e0}, - {-1.18788e-1,-1.33674e0,6.194e0,5.46446e0,-5.62994e1,-7.26915e1}, - {-3.46033e-1,2.03729e-1,-8.38012e0,-8.84466e0,7.34319e1,1.03805e2}, - {6.47255e-2,-1.55409e-1,3.03879e0,3.54516e0,-2.7098e1,-4.01572e1} - }}, - //Br................................................. - CoeffMat{{ - {9.99952e-1,2.79961e-5,4.60479e-4,-8.33486e-4,-2.23214e-3,-5.16285e-4}, - {3.7558e-1,7.68292e-1,-1.62392e0,-4.87674e-1,1.18854e1,1.28025e1}, - {-9.79355e-2,-1.38749e0,7.0556e0,6.80198e0,-6.44113e1,-8.52915e1}, - {-3.66572e-1,2.41863e-1,-9.29524e0,-1.04141e1,8.22093e1,1.18166e2}, - {6.76714e-2,-1.75289e-1,3.33691e0,4.08538e0,-2.99809e1,-4.50381e1} - }}, - //Kr................................................. - CoeffMat{{ - {9.99947e-1,2.2952e-5,6.82639e-4,-7.17139e-4,-4.00522e-3,-2.81601e-3}, - {3.82332e-1,8.08194e-1,-1.86848e0,-8.03415e-1,1.4064e1,1.5956e1}, - {-7.48735e-2,-1.44031e0,7.92233e0,8.2382e0,-7.26858e1,-9.83865e1}, - {-3.88797e-1,2.81837e-1,-1.02005e1,-1.20716e1,9.10175e1,1.32873e2}, - {7.10017e-2,-1.96182e-1,3.6278e0,4.65e0,-3.28365e1,-4.99823e1} - }}, - //Rb................................................. - CoeffMat{{ - {9.99941e-1,1.63607e-5,9.28242e-4,-5.69257e-4,-5.98245e-3,-5.42476e-3}, - {3.88363e-1,8.49548e-1,-2.11701e0,-1.15003e0,1.63119e1,1.92737e1}, - {-4.93364e-2,-1.49491e0,8.78129e0,9.76596e0,-8.10212e1,-1.11851e2}, - {-4.12953e-1,3.23322e-1,-1.10814e1,-1.38073e1,9.97386e1,1.47775e2}, - {7.47911e-2,-2.18001e-1,3.90642e0,5.23509e0,-3.56231e1,-5.49352e1} - }}, - //Sr................................................. - CoeffMat{{ - {9.99935e-1,8.3152e-6,1.19244e-3,-3.95258e-4,-8.12525e-3,-8.28836e-3}, - {3.93603e-1,8.92358e-1,-2.36688e0,-1.52786e0,1.86071e1,2.27314e1}, - {-2.11345e-2,-1.55112e0,9.62205e0,1.13827e1,-8.93266e1,-1.25574e2}, - {-4.39199e-1,3.66161e-1,-1.19262e1,-1.56159e1,1.08267e2,1.62736e2}, - {7.90853e-2,-2.40719e-1,4.16872e0,5.83839e0,-3.8304e1,-5.98484e1} - }}, - //Y................................................. - CoeffMat{{ - {9.99929e-1,-1.65608e-6,1.47476e-3,-1.8625e-4,-1.04363e-2,-1.14318e-2}, - {3.97989e-1,9.3643e-1,-2.61569e0,-1.93387e0,2.09305e1,2.63021e1}, - {9.88744e-3,-1.60813e0,1.04351e1,1.3074e1,-9.75207e1,-1.39437e2}, - {-4.67657e-1,4.09494e-1,-1.2724e1,-1.74797e1,1.16508e2,1.77614e2}, - {8.39169e-2,-2.64077e-1,4.41092e0,6.45344e0,-4.08455e1,-6.46698e1} - }}, - //Zr................................................. - CoeffMat{{ - {9.99922e-1,-1.37624e-5,1.77335e-3,5.994e-5,-1.29013e-2,-1.48443e-2}, - {4.0145e-1,9.816e-1,-2.86053e0,-2.36562e0,2.32591e1,2.99551e1}, - {4.3915e-2,-1.66522e0,1.12093e1,1.48279e1,-1.05512e2,-1.5331e2}, - {-4.98475e-1,4.52607e-1,-1.34628e1,-1.93836e1,1.24357e2,1.92256e2}, - {8.93258e-2,-2.87869e-1,4.62891e0,7.07473e0,-4.32116e1,-6.93461e1} - }}, - //Nb................................................. - CoeffMat{{ - {9.99916e-1,-2.78975e-5,2.08753e-3,3.41607e-4,-1.55141e-2,-1.85156e-2}, - {4.03889e-1,1.02783e0,-3.09762e0,-2.82183e0,2.55629e1,3.36548e1}, - {8.12122e-2,-1.7221e0,1.1931e1,1.66362e1,-1.13185e2,-1.67048e2}, - {-5.3188e-1,4.95236e-1,-1.41275e1,-2.13166e1,1.31687e2,2.06495e2}, - {9.53781e-2,-3.12041e-1,4.81765e0,7.69813e0,-4.53587e1,-7.38182e1} - }}, - //Mo................................................. - CoeffMat{{ - {9.9991e-1,-4.3988e-5,2.41439e-3,6.54636e-4,-1.82489e-2,-2.24062e-2}, - {4.05249e-1,1.07495e0,-3.32445e0,-3.30073e0,2.78221e1,3.7376e1}, - {1.21896e-1,-1.77813e0,1.25907e1,1.84893e1,-1.20461e2,-1.80542e2}, - {-5.67942e-1,5.36741e-1,-1.4708e1,-2.32664e1,1.38408e2,2.20203e2}, - {1.02086e-1,-3.36419e-1,4.97373e0,8.31905e0,-4.72562e1,-7.80409e1} - }}, - //Tc................................................. - CoeffMat{{ - {9.99904e-1,-6.22073e-5,2.75066e-3,1.00063e-3,-2.10817e-2,-2.64933e-2}, - {4.05472e-1,1.12279e0,-3.53825e0,-3.79964e0,3.00132e1,4.10871e1}, - {1.66089e-1,-1.83255e0,1.31785e1,2.03743e1,-1.27252e2,-1.93664e2}, - {-6.06729e-1,5.76392e-1,-1.51936e1,-2.52171e1,1.44423e2,2.33232e2}, - {1.09463e-1,-3.60803e-1,5.09358e0,8.93178e0,-4.88706e1,-8.19632e1} - }}, - //Ru................................................. - CoeffMat{{ - {9.99898e-1,-8.26232e-5,3.09353e-3,1.37924e-3,-2.39902e-2,-3.07502e-2}, - {4.04498e-1,1.17115e0,-3.73628e0,-4.3158e0,3.21142e1,4.4758e1}, - {2.13904e-1,-1.88458e0,1.36848e1,2.22782e1,-1.33471e2,-2.06293e2}, - {-6.48301e-1,6.13439e-1,-1.55742e1,-2.71531e1,1.49638e2,2.45444e2}, - {1.17516e-1,-3.8499e-1,5.17383e0,9.5307e0,-5.01709e1,-8.55368e1} - }}, - //Rh................................................. - CoeffMat{{ - {9.99893e-1,-1.05293e-4,3.43932e-3,1.78981e-3,-2.69462e-2,-3.51436e-2}, - {4.02272e-1,1.21982e0,-3.91575e0,-4.84621e0,3.41013e1,4.83563e1}, - {2.65438e-1,-1.93341e0,1.40998e1,2.41876e1,-1.39034e2,-2.183e2}, - {-6.92691e-1,6.47135e-1,-1.58398e1,-2.90581e1,1.5396e2,2.56697e2}, - {1.26243e-1,-4.08783e-1,5.21125e0,1.011e1,-5.11262e1,-8.8713e1} - }}, - //Pd................................................. - CoeffMat{{ - {9.99888e-1,-1.30188e-4,3.78436e-3,2.22971e-3,-2.99185e-2,-3.96315e-2}, - {3.9875e-1,1.26855e0,-4.07427e0,-5.38796e0,3.5955e1,5.18546e1}, - {3.20748e-1,-1.97817e0,1.44153e1,2.6089e1,-1.43866e2,-2.29578e2}, - {-7.39891e-1,6.76679e-1,-1.59819e1,-3.09162e1,1.57312e2,2.66867e2}, - {1.3563e-1,-4.31972e-1,5.2031e0,1.06642e1,-5.17104e1,-9.14493e1} - }}, - //Ag................................................. - CoeffMat{{ - {9.99884e-1,-1.57459e-4,4.12052e-3,2.69984e-3,-3.28468e-2,-4.41512e-2}, - {3.93891e-1,1.3171e0,-4.20944e0,-5.93691e0,3.76535e1,5.52188e1}, - {3.79857e-1,-2.01791e0,1.46238e1,2.79653e1,-1.47893e2,-2.39997e2}, - {-7.89852e-1,7.01196e-1,-1.59929e1,-3.27077e1,1.59613e2,2.75813e2}, - {1.45646e-1,-4.54328e-1,5.14701e0,1.11863e1,-5.18977e1,-9.36987e1} - }}, - //Cd................................................. - CoeffMat{{ - {9.99881e-1,-1.86816e-4,4.45491e-3,3.19574e-3,-3.57812e-2,-4.87457e-2}, - {3.87612e-1,1.36524e0,-4.31686e0,-6.49005e0,3.91612e1,5.84032e1}, - {4.4294e-1,-2.05189e0,1.47103e1,2.98033e1,-1.50988e2,-2.49389e2}, - {-8.42683e-1,7.20045e-1,-1.58579e1,-3.44169e1,1.60733e2,2.83354e2}, - {1.56312e-1,-4.75705e-1,5.0381e0,1.16708e1,-5.16454e1,-9.53999e1} - }}, - //In................................................. - CoeffMat{{ - {9.99878e-1,-2.1848e-4,4.78376e-3,3.71897e-3,-3.86926e-2,-5.33866e-2}, - {3.79879e-1,1.41265e0,-4.39391e0,-7.04282e0,4.0455e1,6.13721e1}, - {5.10001e-1,-2.07896e0,1.46666e1,3.15843e1,-1.53072e2,-2.57624e2}, - {-8.98306e-1,7.32213e-1,-1.55686e1,-3.60229e1,1.60591e2,2.89348e2}, - {1.67591e-1,-4.95837e-1,4.87389e0,1.21106e1,-5.09274e1,-9.6506e1} - }}, - //Sn................................................ - CoeffMat{{ - {9.99876e-1,-2.52173e-4,5.09558e-3,4.26518e-3,-4.14939e-2,-5.79712e-2}, - {3.70689e-1,1.45908e0,-4.43983e0,-7.59192e0,4.15261e1,6.4108e1}, - {5.80924e-1,-2.09831e0,1.44911e1,3.32945e1,-1.54118e2,-2.64634e2}, - {-9.56518e-1,7.37009e-1,-1.51243e1,-3.75098e1,1.59159e2,2.93723e2}, - {1.79397e-1,-5.14573e-1,4.65426e0,1.25003e1,-4.97362e1,-9.69932e1} - }}, - //Sb................................................. - CoeffMat{{ - {9.99874e-1,-2.87917e-4,5.39066e-3,4.83313e-3,-4.41852e-2,-6.24953e-2}, - {3.60001e-1,1.5042e0,-4.4515e0,-8.1327e0,4.23476e1,6.65709e1}, - {6.55716e-1,-2.10882e0,1.41741e1,3.49157e1,-1.54034e2,-2.70276e2}, - {-1.01724e0,7.33501e-1,-1.45156e1,-3.88572e1,1.56348e2,2.96328e2}, - {1.9169e-1,-5.3169e-1,4.37639e0,1.28328e1,-4.80435e1,-9.68125e1} - }}, - //Te................................................. - CoeffMat{{ - {9.99874e-1,-3.25847e-4,5.67083e-3,5.42945e-3,-4.67877e-2,-6.701e-2}, - {3.47783e-1,1.54767e0,-4.42606e0,-8.65995e0,4.28934e1,6.87183e1}, - {7.34353e-1,-2.10941e0,1.3707e1,3.64273e1,-1.52736e2,-2.74401e2}, - {-1.08036e0,7.20771e-1,-1.37344e1,-4.00427e1,1.52073e2,2.97007e2}, - {2.04417e-1,-5.46972e-1,4.03784e0,1.31007e1,-4.58224e1,-9.59123e1} - }}, - //I................................................. - CoeffMat{{ - {9.99875e-1,-3.65406e-4,5.92115e-3,6.03332e-3,-4.91674e-2,-7.12951e-2}, - {3.34062e-1,1.58924e0,-4.36361e0,-9.17241e0,4.31638e1,7.05522e1}, - {8.16604e-1,-2.09933e0,1.30915e1,3.78235e1,-1.50231e2,-2.77015e2}, - {-1.14554e0,6.98298e-1,-1.2784e1,-4.10591e1,1.46348e2,2.9577e2}, - {2.17451e-1,-5.60344e-1,3.63994e0,1.33015e1,-4.30802e1,-9.42981e1} - }}, - //Xe.................................................. - CoeffMat{{ - {9.99877e-1,-4.06637e-4,6.14278e-3,6.64937e-3,-5.13391e-2,-7.53881e-2}, - {3.18822e-1,1.62854e0,-4.26169e0,-9.66459e0,4.31357e1,7.20331e1}, - {9.02373e-1,-2.07743e0,1.23209e1,3.90832e1,-1.46445e2,-2.77984e2}, - {-1.21259e0,6.65174e-1,-1.16583e1,-4.18839e1,1.39104e2,2.92479e2}, - {2.3071e-1,-5.71608e-1,3.18104e0,1.34277e1,-3.97963e1,-9.19255e1} - }}, - //Cs.................................................. - CoeffMat{{ - {9.99881e-1,-4.49175e-4,6.32631e-3,7.27184e-3,-5.32297e-2,-7.91972e-2}, - {3.02086e-1,1.66529e0,-4.11996e0,-1.0133e1,4.28031e1,7.31464e1}, - {9.91429e-1,-2.04292e0,1.13959e1,4.01931e1,-1.41369e2,-2.7726e2}, - {-1.28117e0,6.2087e-1,-1.03596e1,-4.25024e1,1.30338e2,2.8709e2}, - {2.44065e-1,-5.80702e-1,2.66217e0,1.34745e1,-3.59721e1,-8.87824e1} - }}, - //Ba.................................................. - CoeffMat{{ - {9.99886e-1,-4.93136e-4,6.47914e-3,7.90029e-3,-5.48909e-2,-8.27736e-2}, - {2.8384e-1,1.69905e0,-3.93534e0,-1.05724e1,4.21392e1,7.38521e1}, - {1.08366e0,-1.99451e0,1.03076e1,4.1134e1,-1.34917e2,-2.7471e2}, - {-1.35107e0,5.64397e-1,-8.87993e0,-4.28945e1,1.19972e2,2.7947e2}, - {2.57431e-1,-5.87417e-1,2.08111e0,1.34353e1,-3.15844e1,-8.48268e1} - }}, - //La.................................................. - CoeffMat{{ - {9.99892e-1,-5.38122e-4,6.59206e-3,8.52848e-3,-5.62508e-2,-8.60251e-2}, - {2.6413e-1,1.72952e0,-3.70812e0,-1.09794e1,4.1143e1,7.4141e1}, - {1.17876e0,-1.93141e0,9.05928e0,4.18928e1,-1.27099e2,-2.70309e2}, - {-1.42185e0,4.95292e-1,-7.22427e0,-4.30461e1,1.08024e2,2.69601e2}, - {2.70647e-1,-5.91727e-1,1.43982e0,1.33056e1,-2.66423e1,-8.00557e1} - }}, - //Ce.................................................. - CoeffMat{{ - {9.999e-1,-5.83896e-4,6.66373e-3,9.15389e-3,-5.72988e-2,-8.89368e-2}, - {2.42961e-1,1.75635e0,-3.43619e0,-1.13482e1,3.97928e1,7.39747e1}, - {1.27654e0,-1.85263e0,7.64589e0,4.24492e1,-1.17853e2,-2.63933e2}, - {-1.49322e0,4.12926e-1,-5.38909e0,-4.29363e1,9.44419e1,2.57363e2}, - {2.83597e-1,-5.93557e-1,7.3754e-1,1.30785e1,-2.11316e1,-7.44329e1} - }}, - //Pr.................................................. - CoeffMat{{ - {9.99909e-1,-6.301e-4,6.69295e-3,9.77077e-3,-5.8022e-2,-9.14799e-2}, - {2.20373e-1,1.77919e0,-3.11907e0,-1.16754e1,3.80818e1,7.33384e1}, - {1.37668e0,-1.75734e0,6.06819e0,4.27905e1,-1.07169e2,-2.55543e2}, - {-1.56475e0,3.16885e-1,-3.37687e0,-4.25518e1,7.92269e1,2.42726e2}, - {2.96124e-1,-5.92915e-1,-2.45072e-2,1.27501e1,-1.50562e1,-6.79514e1} - }}, - //Nd.................................................. - CoeffMat{{ - {9.9992e-1,-6.77211e-4,6.67917e-3,1.03866e-2,-5.84208e-2,-9.36775e-2}, - {1.96406e-1,1.79753e0,-2.7553e0,-1.19533e1,3.59929e1,7.21943e1}, - {1.47887e0,-1.64416e0,4.32381e0,4.28906e1,-9.50034e1,-2.45021e2}, - {-1.636e0,2.06169e-1,-1.18704e0,-4.18661e1,6.23476e1,2.25579e2}, - {3.08067e-1,-5.89622e-1,-8.4608e-1,1.23116e1,-8.40939e0,-6.05788e1} - }}, - //Pm.................................................. - CoeffMat{{ - {9.99932e-1,-7.23543e-4,6.61026e-3,1.09756e-2,-5.83892e-2,-9.53462e-2}, - {1.71135e-1,1.81127e0,-2.34665e0,-1.21818e1,3.35368e1,7.05547e1}, - {1.58269e0,-1.5131e0,2.4214e0,4.27487e1,-8.14091e1,-2.32422e2}, - {-1.70642e0,8.13302e-2,1.16945e0,-4.0879e1,4.38728e1,2.05995e2}, - {3.19226e-1,-5.84023e-1,-1.72313e0,1.17636e1,-1.21776e0,-5.23434e1} - }}, - //Sm.................................................. - CoeffMat{{ - {9.99946e-1,-7.69368e-4,6.48923e-3,1.15425e-2,-5.79523e-2,-9.65284e-2}, - {1.44613e-1,1.81998e0,-1.89196e0,-1.23548e1,3.06998e1,6.83894e1}, - {1.68776e0,-1.36316e0,3.59801e-1,4.23443e1,-6.63559e1,-2.17658e2}, - {-1.77552e0,-5.81387e-2,3.69184e0,-3.95701e1,2.37855e1,1.83897e2}, - {3.2942e-1,-5.76119e-1,-2.65496e0,1.10995e1,6.52027e0,-4.32241e1} - }}, - //Eu.................................................. - CoeffMat{{ - {9.99962e-1,-8.13654e-4,6.31217e-3,1.20763e-2,-5.70782e-2,-9.71683e-2}, - {1.16915e-1,1.82343e0,-1.39208e0,-1.24689e1,2.74833e1,6.56914e1}, - {1.79364e0,-1.194e0,-1.85513e0,4.16665e1,-4.9868e1,-2.00724e2}, - {-1.84272e0,-2.12012e-1,6.37189e0,-3.79293e1,2.12742e0,1.59295e2}, - {3.38449e-1,-5.66167e-1,-3.63838e0,1.03165e1,1.47866e1,-3.32296e1} - }}, - //Gd.................................................. - CoeffMat{{ - {9.9998e-1,-8.57331e-4,6.07736e-3,1.25819e-2,-5.57471e-2,-9.72475e-2}, - {8.8154e-2,1.82105e0,-8.47718e-1,-1.25186e1,2.38891e1,6.2451e1}, - {1.89976e0,-1.00428e0,-4.21799e0,4.0697e1,-3.19703e1,-1.81604e2}, - {-1.90731e0,-3.8108e-1,9.2019e0,-3.59386e1,-2.10594e1,1.32192e2}, - {3.46059e-1,-5.54093e-1,-4.67039e0,9.40924e0,2.35631e1,-2.23653e1} - }}, - //Tb.................................................. - CoeffMat{{ - {9.99999e-1,-8.99343e-4,5.7868e-3,1.30512e-2,-5.39764e-2,-9.67723e-2}, - {5.84024e-2,1.81266e0,-2.59193e-1,-1.25002e1,1.99139e1,5.86545e1}, - {2.00567e0,-7.93885e-1,-6.72461e0,3.94244e1,-1.26714e1,-1.60273e2}, - {-1.96873e0,-5.64782e-1,1.21753e1,-3.3588e1,-4.57487e1,1.02578e2}, - {3.5205e-1,-5.40282e-1,-5.74837e0,8.37487e0,3.28364e1,-1.06341e1} - }}, - //Dy.................................................. - CoeffMat{{ - {1.00002e0,-9.39311e-4,5.43607e-3,1.34763e-2,-5.17249e-2,-9.56765e-2}, - {2.78005e-2,1.79795e0,3.70737e-1,-1.24114e1,1.55759e1,5.43184e1}, - {2.1107e0,-5.62473e-1,-9.36239e0,3.78424e1,7.94516e0,-1.36806e2}, - {-2.02617e0,-7.62769e-1,1.52768e1,-3.08728e1,-7.18354e1,7.05554e1}, - {3.5614e-1,-5.25061e-1,-6.86679e0,7.21247e0,4.25672e1,1.92463e0} - }}, - //Ho................................................... - CoeffMat{{ - {1.00004e0,-9.76747e-4,5.02266e-3,1.3853e-2,-4.89737e-2,-9.3934e-2}, - {-3.54744e-3,1.77664e0,1.04132e0,-1.22473e1,1.08741e1,4.94269e1}, - {2.21429e0,-3.09743e-1,-1.21254e1,3.59359e1,2.986e1,-1.11177e2}, - {-2.07893e0,-9.74603e-1,1.84979e1,-2.77795e1,-9.928e1,3.61173e1}, - {3.5809e-1,-5.08796e-1,-8.02233e0,5.91833e0,5.27374e1,1.53069e1} - }}, - //Er................................................... - CoeffMat{{ - {1.00007e0,-1.01117e-3,4.54349e-3,1.4172e-2,-4.56907e-2,-9.14834e-2}, - {-3.54805e-2,1.74848e0,1.74921e0,-1.20063e1,5.83106e0,4.40029e1}, - {2.3157e0,-3.56472e-2,-1.49988e1,3.37013e1,5.29727e1,-8.34879e1}, - {-2.12613e0,-1.19953e0,2.18208e1,-2.43067e1,-1.27958e2,-6.05797e-1}, - {3.57593e-1,-4.91966e-1,-9.2087e0,4.49274e0,6.33009e1,2.94631e1} - }}, - //Tm................................................... - CoeffMat{{ - {1.00009e0,-1.04219e-3,3.99979e-3,1.44301e-2,-4.18849e-2,-8.83308e-2}, - {-6.78832e-2,1.71323e0,2.49363e0,-1.16834e1,4.45978e-1,3.80313e1}, - {2.41433e0,2.59886e-1,-1.79765e1,3.11249e1,7.72617e1,-5.37152e1}, - {-2.16704e0,-1.4368e0,2.52369e1,-2.04429e1,-1.57828e2,-3.9614e1}, - {3.544e-1,-4.75067e-1,-1.04225e1,2.93269e0,7.42385e1,4.43866e1} - }}, - //Yb................................................... - CoeffMat{{ - {1.00012e0,-1.06927e-3,3.39091e-3,1.46208e-2,-3.75471e-2,-8.44506e-2}, - {-1.00611e-1,1.67068e0,3.27246e0,-1.12756e1,-5.27058e0,3.15152e1}, - {2.50948e0,5.76648e-1,-2.10477e1,2.81988e1,1.02665e2,-2.18989e1}, - {-2.20084e0,-1.68528e0,2.87325e1,-1.61828e1,-1.88802e2,-8.08394e1}, - {3.48226e-1,-4.58722e-1,-1.16589e1,1.23735e0,8.55166e1,6.00477e1} - }}, - //Lu................................................... - CoeffMat{{ - {1.00015e0,-1.09216e-3,2.71745e-3,1.47404e-2,-3.26779e-2,-7.98343e-2}, - {-1.33495e-1,1.62059e0,4.08307e0,-1.07793e1,-1.13042e1,2.44607e1}, - {2.60038e0,9.14519e-1,-2.42e1,2.49142e1,1.29108e2,1.19101e1}, - {-2.22661e0,-1.94391e0,3.22924e1,-1.15206e1,-2.20784e2,-1.24202e2}, - {3.38764e-1,-4.43546e-1,-1.29122e1,-5.9428e-1,9.70972e1,7.64128e1} - }}, - //Hf.................................................. - CoeffMat{{ - {1.00018e0,-1.11005e-3,1.98019e-3,1.47814e-2,-2.72814e-2,-7.44728e-2}, - {-1.66392e-1,1.56288e0,4.92359e0,-1.01915e1,-1.76472e1,1.68656e1}, - {2.68631e0,1.27279e0,-2.74231e1,2.12642e1,1.56537e2,4.76851e1}, - {-2.24352e0,-2.21096e0,3.59035e1,-6.45262e0,-2.53694e2,-1.69646e2}, - {3.25734e-1,-4.30391e-1,-1.41779e1,-2.56236e0,1.08949e2,9.34555e1} - }}, - //Ta.................................................. - CoeffMat{{ - {1.00021e0,-1.1227e-3,1.1875e-3,1.47467e-2,-2.14258e-2,-6.84575e-2}, - {-1.99145e-1,1.49739e0,5.79226e0,-9.50812e0,-2.42938e1,8.72312e0}, - {2.76656e0,1.65088e0,-3.07077e1,1.72391e1,1.84905e2,8.54178e1}, - {-2.25073e0,-2.48478e0,3.95541e1,-9.72452e-1,-2.87462e2,-2.17133e2}, - {3.0885e-1,-4.20088e-1,-1.54515e1,-4.66797e0,1.21043e2,1.11156e2} - }}, - //W................................................... - CoeffMat{{ - {1.00025e0,-1.12974e-3,3.25083e-4,1.46196e-2,-1.49818e-2,-6.15916e-2}, - {-2.31499e-1,1.42403e0,6.68301e0,-8.72894e0,-3.12005e1,8.13039e-2}, - {2.84004e0,2.04799e0,-3.40292e1,1.28426e1,2.14035e2,1.24907e2}, - {-2.24701e0,-2.76343e0,4.3216e1,4.91177e0,-3.2188e2,-2.66423e2}, - {2.87708e-1,-4.13579e-1,-1.67235e1,-6.90724e0,1.33305e2,1.29426e2} - }}, - //Re.................................................. - CoeffMat{{ - {1.00028e0,-1.13018e-3,-5.85208e-4,1.44008e-2,-8.12721e-3,-5.40793e-2}, - {-2.63328e-1,1.34285e0,7.59541e0,-7.85118e0,-3.83723e1,-9.07534e0}, - {2.90612e0,2.46267e0,-3.73829e1,8.06978e0,2.43918e2,1.66172e2}, - {-2.23164e0,-3.04426e0,4.68824e1,1.12007e1,-3.56918e2,-3.17506e2}, - {2.62068e-1,-4.12056e-1,-1.79911e1,-9.27922e0,1.45719e2,1.48255e2} - }}, - //Os.................................................. - CoeffMat{{ - {1.00032e0,-1.12442e-3,-1.55891e-3,1.40881e-2,-7.3226e-4,-4.57754e-2}, - {-2.94323e-1,1.25374e0,8.52167e0,-6.87354e0,-4.57537e1,-1.8691e1}, - {2.96357e0,2.89403e0,-4.07385e1,2.92128e0,2.74335e2,2.08986e2}, - {-2.20324e0,-3.32515e0,5.05197e1,1.78893e1,-3.92323e2,-3.70113e2}, - {2.31474e-1,-4.16536e-1,-1.92429e1,-1.1781e1,1.58197e2,1.67545e2} - }}, - //Ir................................................... - CoeffMat{{ - {1.00036e0,-1.11107e-3,-2.58076e-3,1.36692e-2,7.08911e-3,-3.6773e-2}, - {-3.24316e-1,1.15695e0,9.45947e0,-5.79675e0,-5.33328e1,-2.87528e1}, - {3.01159e0,3.33981e0,-4.4085e1,-2.59468e0,3.05217e2,2.53266e2}, - {-2.16093e0,-3.60244e0,5.41141e1,2.49638e1,-4.28002e2,-4.24126e2}, - {1.95645e-1,-4.28567e-1,-2.04742e1,-1.44066e1,1.70702e2,1.87249e2} - }}, - //Pt................................................... - CoeffMat{{ - {1.00039e0,-1.09011e-3,-3.65347e-3,1.31413e-2,1.53649e-2,-2.70335e-2}, - {-3.53028e-1,1.05257e0,1.04028e1,-4.62034e0,-6.10681e1,-3.92186e1}, - {3.04907e0,3.79831e0,-4.73979e1,-8.47391e0,3.36394e2,2.98829e2}, - {-2.10347e0,-3.87306e0,5.76385e1,3.24147e1,-4.63751e2,-4.79318e2}, - {1.54174e-1,-4.49506e-1,-2.16756e1,-1.71516e1,1.83161e2,2.07283e2} - }}, - //Au................................................... - CoeffMat{{ - {1.00043e0,-1.06114e-3,-4.7734e-3,1.25001e-2,2.40713e-2,-1.65728e-2}, - {-3.80207e-1,9.40778e-1,1.13464e1,-3.34423e0,-6.89259e1,-5.00549e1}, - {3.07493e0,4.26732e0,-5.06567e1,-1.47102e1,3.67721e2,3.45521e2}, - {-2.0297e0,-4.13334e0,6.10696e1,4.02298e1,-4.99399e2,-5.35495e2}, - {1.06694e-1,-4.80932e-1,-2.28393e1,-2.00106e1,1.95512e2,2.27574e2} - }}, - //Hg................................................. - CoeffMat{{ - {1.00047e0,-1.02377e-3,-5.93877e-3,1.17393e-2,3.32021e-2,-5.38025e-3}, - {-4.05577e-1,8.21834e-1,1.22846e1,-1.96909e0,-7.68675e1,-6.12207e1}, - {3.08808e0,4.74434e0,-5.38383e1,-2.12945e1,3.99035e2,3.93158e2}, - {-1.93842e0,-4.37923e0,6.43816e1,4.83938e1,-5.3475e2,-5.92428e2}, - {5.28167e-2,-5.24544e-1,-2.39565e1,-2.29768e1,2.07685e2,2.48037e2} - }}, - //Tl................................................. - CoeffMat{{ - {1.00051e0,-9.77691e-4,-7.14416e-3,1.08563e-2,4.27196e-2,6.50699e-3}, - {-4.28882e-1,6.96049e-1,1.32126e1,-4.94649e-1,-8.48638e1,-7.26907e1}, - {3.08745e0,5.22659e0,-5.69235e1,-2.82201e1,4.30208e2,4.41611e2}, - {-1.82849e0,-4.60634e0,6.7553e1,5.68935e1,-5.69648e2,-6.49946e2}, - {-7.81506e-3,-5.82177e-1,-2.50201e1,-2.60444e1,2.19621e2,2.68607e2} - }}, - //Pb................................................. - CoeffMat{{ - {1.00055e0,-9.22411e-4,-8.38967e-3,9.83942e-3,5.26394e-2,1.91419e-2}, - {-4.49833e-1,5.63865e-1,1.41246e1,1.0768e0,-9.28737e1,-8.44177e1}, - {3.07189e0,5.71077e0,-5.98888e1,-3.54715e1,4.61066e2,4.9067e2}, - {-1.69869e0,-4.80968e0,7.05579e1,6.57062e1,-6.03887e2,-7.07793e2}, - {-7.55851e-2,-6.55871e-1,-2.60213e1,-2.92042e1,2.31249e2,2.89188e2} - }}, - //Bi................................................. - CoeffMat{{ - {1.00059e0,-8.58152e-4,-9.6602e-3,8.70229e-3,6.28328e-2,3.23376e-2}, - {-4.68192e-1,4.25668e-1,1.50174e1,2.74763e0,-1.00883e2,-9.64012e1}, - {3.04043e0,6.19366e0,-6.27207e1,-4.30487e1,4.91532e2,5.40291e2}, - {-1.54799e0,-4.98431e0,7.3381e1,7.48252e1,-6.37365e2,-7.65885e2}, - {-1.5081e-1,-7.47659e-1,-2.69551e1,-3.24524e1,2.42528e2,3.09744e2} - }}, - //Po................................................. - CoeffMat{{ - {1.00062e0,-7.84582e-4,-1.09553e-2,7.43533e-3,7.33106e-2,4.61358e-2}, - {-4.83646e-1,2.8205e-1,1.58842e1,4.51493e0,-1.08846e2,-1.08589e2}, - {2.99185e0,6.67128e0,-6.53938e1,-5.09329e1,5.21419e2,5.90246e2}, - {-1.37512e0,-5.12439e0,7.59949e1,8.42238e1,-6.69862e2,-8.23945e2}, - {-2.33878e-1,-8.59879e-1,-2.78125e1,-3.57783e1,2.53381e2,3.30173e2} - }}, - //At................................................. - CoeffMat{{ - {1.00066e0,-7.01484e-4,-1.22684e-2,6.03454e-3,8.40298e-2,6.05018e-2}, - {-4.95899e-1,1.3369e-1,1.67195e1,6.37574e0,-1.16723e2,-1.20936e2}, - {2.92499e0,7.13921e0,-6.78861e1,-5.91051e1,5.50562e2,6.40337e2}, - {-1.17889e0,-5.22361e0,7.8376e1,9.38748e1,-7.01185e2,-8.81728e2}, - {-3.25146e-1,-9.95032e-1,-2.85858e1,-3.91713e1,2.63738e2,3.50385e2} - }}, - //Rn................................................. - CoeffMat{{ - {1.0007e0,-6.08729e-4,-1.35873e-2,4.50061e-3,9.48991e-2,7.53277e-2}, - {-5.04647e-1,-1.86445e-2,1.75175e1,8.32729e0,-1.24477e2,-1.33405e2}, - {2.83871e0,7.59269e0,-7.01766e1,-6.75463e1,5.78809e2,6.90388e2}, - {-9.58143e-1,-5.27523e0,8.0502e1,1.03751e2,-7.31155e2,-9.39013e2}, - {-4.24964e-1,-1.15578e0,-2.92679e1,-4.26205e1,2.73537e2,3.70298e2} - }}, - //Fr................................................. - CoeffMat{{ - {1.00073e0,-5.06594e-4,-1.49195e-2,2.82752e-3,1.05984e-1,9.07179e-2}, - {-5.09498e-1,-1.74114e-1,1.82691e1,1.0364e1,-1.32042e2,-1.45912e2}, - {2.73154e0,8.02658e0,-7.22315e1,-7.62278e1,6.05902e2,7.40063e2}, - {-7.11399e-1,-5.27209e0,8.23375e1,1.13815e2,-7.59479e2,-9.95409e2}, - {-5.33773e-1,-1.34491e0,-2.98474e1,-4.61115e1,2.82675e2,3.89772e2} - }}, - //Ra................................................. - CoeffMat{{ - {1.00076e0,-3.94938e-4,-1.62461e-2,1.01616e-3,1.17156e-1,1.06529e-1}, - {-5.10179e-1,-3.3169e-1,1.89693e1,1.24811e1,-1.39382e2,-1.58412e2}, - {2.60247e0,8.43503e0,-7.4033e1,-8.51228e1,6.31692e2,7.89166e2}, - {-4.37673e-1,-5.2062e0,8.38642e1,1.2403e2,-7.85984e2,-1.05067e3}, - {-6.5185e-1,-1.56551e0,-3.03186e1,-4.96307e1,2.91091e2,4.08719e2} - }}, - //Ac................................................. - CoeffMat{{ - {1.00079e0,-2.74014e-4,-1.75614e-2,-9.33046e-4,1.28369e-1,1.22714e-1}, - {-5.06317e-1,-4.90282e-1,1.96102e1,1.46725e1,-1.4644e2,-1.70838e2}, - {2.45014e0,8.81187e0,-7.5553e1,-9.42005e1,6.55963e2,8.37419e2}, - {-1.35635e-1,-5.06921e0,8.50532e1,1.34356e2,-8.10425e2,-1.10448e3}, - {-7.79579e-1,-1.82077e0,-3.06725e1,-5.31628e1,2.987e2,4.2702e2} - }}, - //Th................................................ - CoeffMat{{ - {1.00082e0,-1.43895e-4,-1.8859e-2,-3.02404e-3,1.39581e-1,1.39236e-1}, - {-4.97537e-1,-6.48585e-1,2.01838e1,1.69296e1,-1.53157e2,-1.83108e2}, - {2.27321e0,9.15011e0,-7.6763e1,-1.0342e2,6.78484e2,8.84496e2}, - {1.9599e-1,-4.85189e0,8.5876e1,1.44742e2,-8.32543e2,-1.15645e3}, - {-9.17319e-1,-2.1142e0,-3.09003e1,-5.66889e1,3.05413e2,4.44542e2} - }}, - //Pa................................................. - CoeffMat{{ - {1.00085e0,-5.30411e-6,-2.01239e-2,-5.2478e-3,1.50683e-1,1.55963e-1}, - {-4.83526e-1,-8.05262e-1,2.06846e1,1.92461e1,-1.59496e2,-1.95177e2}, - {2.07059e0,9.44256e0,-7.76442e1,-1.12749e2,6.99103e2,9.30193e2}, - {5.58224e-1,-4.54471e0,8.63141e1,1.55146e2,-8.52163e2,-1.20633e3}, - {-1.06534e0,-2.44938e0,-3.09967e1,-6.01932e1,3.11168e2,4.61192e2} - }}, - //U.................................................... - CoeffMat{{ - {1.00087e0,1.42049e-4,-2.13455e-2,-7.61031e-3,1.61591e-1,1.72814e-1}, - {-4.6392e-1,-9.58725e-1,2.11056e1,2.16134e1,-1.65402e2,-2.06974e2}, - {1.841e0,9.68106e0,-7.81727e1,-1.22145e2,7.17614e2,9.74218e2}, - {9.52222e-1,-4.13715e0,8.63446e1,1.65515e2,-8.69054e2,-1.25378e3}, - {-1.22395e0,-2.83024e0,-3.09548e1,-6.36564e1,3.15887e2,4.76849e2} - }}}; - // clang-format on - static_assert(std::size(mott_coeffs) - == CoulombScatteringElementData::num_mott_elements, - "wrong number of Mott coefficient elements"); - - int index = z.unchecked_get() - 1; - CELER_VALIDATE( - index >= 0 - && index < int{CoulombScatteringElementData::num_mott_elements}, - << "atomic number " << z.get() - << " is out of range for Coulomb scattering model Mott coefficients " - "(must be less than " - << CoulombScatteringElementData::num_mott_elements << ")"); - - return mott_coeffs[index]; -} - -//---------------------------------------------------------------------------// -/*! - * Calculate the constant prefactors of the squared momentum transfer. - * - * This factor is used in the exponential and Gaussian nuclear form models: see - * Eqs. 2.262--2.264 of [LR15]. - * - * Specifically, it calculates \f$ (r_n/\bar h)^2 / 12 \f$. A special case is - * inherited from Geant for hydrogen targets. - */ -real_type -CoulombScatteringModel::calc_nuclear_form_prefactor(IsotopeView const& iso) -{ - if (iso.atomic_number().get() == 1) - { - // TODO: Geant4 hardcodes a different prefactor for hydrogen - return real_type{1.5485e-6}; - } - - // The ratio has units of (MeV/c)^-2, so it's easier to convert the - // inverse which has units of MomentumSq, then invert afterwards - constexpr real_type ratio - = 1 - / native_value_to( - 12 - * ipow<2>(constants::hbar_planck - / (real_type(1.27) * units::femtometer))) - .value(); - return ratio - * fastpow(real_type(iso.atomic_mass_number().get()), - 2 * real_type(0.27)); + return this->host_ref().action; } //---------------------------------------------------------------------------// diff --git a/src/celeritas/em/model/CoulombScatteringModel.cu b/src/celeritas/em/model/CoulombScatteringModel.cu index 6cb08ffc1f..d6614b925a 100644 --- a/src/celeritas/em/model/CoulombScatteringModel.cu +++ b/src/celeritas/em/model/CoulombScatteringModel.cu @@ -8,6 +8,7 @@ #include "CoulombScatteringModel.hh" #include "celeritas/em/executor/CoulombScatteringExecutor.hh" +#include "celeritas/em/params/WentzelOKVIParams.hh" #include "celeritas/global/ActionLauncher.device.hh" #include "celeritas/global/CoreParams.hh" #include "celeritas/global/CoreState.hh" @@ -23,11 +24,14 @@ namespace celeritas void CoulombScatteringModel::execute(CoreParams const& params, CoreStateDevice& state) const { + CELER_EXPECT(params.wentzel()); + auto execute = make_action_track_executor( params.ptr(), state.ptr(), this->action_id(), - InteractionApplier{CoulombScatteringExecutor{this->device_ref()}}); + InteractionApplier{CoulombScatteringExecutor{ + this->device_ref(), params.wentzel()->device_ref()}}); static ActionLauncher const launch_kernel(*this); launch_kernel(state, execute); } diff --git a/src/celeritas/em/model/CoulombScatteringModel.hh b/src/celeritas/em/model/CoulombScatteringModel.hh index 236001c05b..af56db6706 100644 --- a/src/celeritas/em/model/CoulombScatteringModel.hh +++ b/src/celeritas/em/model/CoulombScatteringModel.hh @@ -9,11 +9,8 @@ #include -#include "corecel/data/CollectionMirror.hh" #include "celeritas/em/data/CoulombScatteringData.hh" -#include "celeritas/phys/AtomicNumber.hh" #include "celeritas/phys/ImportedModelAdapter.hh" -#include "celeritas/phys/ImportedProcessAdapter.hh" #include "celeritas/phys/Model.hh" namespace celeritas @@ -32,33 +29,12 @@ class CoulombScatteringModel final : public Model //!@{ //! \name Type aliases using SPConstImported = std::shared_ptr; - using HostRef = CoulombScatteringHostRef; - using DeviceRef = CoulombScatteringDeviceRef; //!@} - //! Wentzel Coulomb scattering model configuration options - struct Options - { - //! Nuclear form factor model - NuclearFormFactorType form_factor_model{ - NuclearFormFactorType::exponential}; - - //! User defined screening factor - real_type screening_factor{1}; - - //! Whether to use integral method to sample interaction length - bool use_integral_xs{true}; - - //! Check if the options are valid - explicit operator bool() const { return screening_factor > 0; } - }; - public: // Construct from model ID and other necessary data CoulombScatteringModel(ActionId id, ParticleParams const& particles, - MaterialParams const& materials, - Options const& options, SPConstImported data); // Particle types and energy ranges that this model applies to @@ -87,24 +63,13 @@ class CoulombScatteringModel final : public Model //!@{ //! Access model data - HostRef const& host_ref() const { return data_.host_ref(); } - DeviceRef const& device_ref() const { return data_.device_ref(); } + CoulombScatteringData const& host_ref() const { return data_; } + CoulombScatteringData const& device_ref() const { return data_; } //!@} private: - CollectionMirror data_; + CoulombScatteringData data_; ImportedModelAdapter imported_; - - // Construct per element data (loads Mott coefficients) - void build_data(HostVal& host_data, - MaterialParams const& materials); - - // Retrieve matrix of interpolated Mott coefficients - static CoulombScatteringElementData::MottCoeffMatrix - get_mott_coeff_matrix(AtomicNumber z); - - // Calculate the nuclear form prefactor - static real_type calc_nuclear_form_prefactor(IsotopeView const& iso); }; //---------------------------------------------------------------------------// diff --git a/src/celeritas/em/params/UrbanMscParams.cc b/src/celeritas/em/params/UrbanMscParams.cc index 83256eb37d..214318ae9b 100644 --- a/src/celeritas/em/params/UrbanMscParams.cc +++ b/src/celeritas/em/params/UrbanMscParams.cc @@ -43,10 +43,7 @@ UrbanMscParams::from_import(ParticleParams const& particles, MaterialParams const& materials, ImportData const& data) { - auto is_urban = [](ImportMscModel const& imm) { - return imm.model_class == ImportModelClass::urban_msc; - }; - if (!std::any_of(data.msc_models.begin(), data.msc_models.end(), is_urban)) + if (!has_msc_model(data, ImportModelClass::urban_msc)) { // No Urban MSC present return nullptr; diff --git a/src/celeritas/em/params/WentzelOKVIParams.cc b/src/celeritas/em/params/WentzelOKVIParams.cc new file mode 100644 index 0000000000..f21043f116 --- /dev/null +++ b/src/celeritas/em/params/WentzelOKVIParams.cc @@ -0,0 +1,951 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 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/params/WentzelOKVIParams.cc +//---------------------------------------------------------------------------// +#include "WentzelOKVIParams.hh" + +#include +#include + +#include "corecel/io/Logger.hh" +#include "corecel/sys/ScopedMem.hh" +#include "celeritas/io/ImportData.hh" +#include "celeritas/mat/MaterialParams.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Construct if Wentzel VI or Coulomb is present, else return nullptr. + */ +std::shared_ptr +WentzelOKVIParams::from_import(ImportData const& data, + SPConstMaterials materials) +{ + CELER_EXPECT(materials); + + bool wentzel = has_msc_model(data, ImportModelClass::wentzel_vi_uni); + bool coulomb = has_model(data, ImportModelClass::e_coulomb_scattering); + if (!(wentzel || coulomb)) + { + // No Wentzel VI MSC or Coulomb scattering present + return nullptr; + } + + Options opts; + opts.is_combined = wentzel && coulomb; + opts.polar_angle_limit = [&]() -> real_type { + if (!coulomb) + { + // Set the maximum scattering angle for Wentzel VI MSC + return constants::pi; + } + if (!wentzel) + { + // Set the minimum scattering angle for Coulomb single scattering + return real_type(0); + } + // Polar angle limit between single and multiple scattering if both + // models are present + return data.em_params.msc_theta_limit; + }(); + opts.screening_factor = data.em_params.screening_factor; + opts.angle_limit_factor = data.em_params.angle_limit_factor; + opts.form_factor = data.em_params.form_factor; + + return std::make_shared(materials, opts); +} + +//---------------------------------------------------------------------------// +/*! + * Construct from cross section data and material properties. + */ +WentzelOKVIParams::WentzelOKVIParams(SPConstMaterials materials, + Options options) +{ + CELER_EXPECT(materials); + + ScopedMem record_mem("WentzelOKVIParams.construct"); + + HostVal host_data; + + host_data.params.is_combined = options.is_combined; + host_data.params.costheta_limit = std::cos(options.polar_angle_limit); + host_data.params.a_sq_factor + = real_type(0.5) + * ipow<2>(options.angle_limit_factor * constants::hbar_planck + * constants::c_light * units::femtometer); + host_data.params.screening_factor = options.screening_factor; + host_data.params.form_factor_type = options.form_factor; + + // Load Mott coefficients + build_data(host_data, *materials); + CELER_ASSERT(host_data); + + // Move to mirrored data, copying to device + data_ = CollectionMirror{std::move(host_data)}; + CELER_ENSURE(data_); +} + +//---------------------------------------------------------------------------// +/*! + * Load Mott coefficients and construct per-element data. + */ +void WentzelOKVIParams::build_data(HostVal& host_data, + MaterialParams const& materials) +{ + // Build element data + size_type const num_elements = materials.num_elements(); + auto elem_data = make_builder(&host_data.elem_data); + elem_data.reserve(num_elements); + + for (auto el_id : range(ElementId{num_elements})) + { + // Load Mott coefficients + MottElementData z_data; + z_data.mott_coeff + = get_mott_coeff_matrix(materials.get(el_id).atomic_number()); + elem_data.push_back(z_data); + } + + auto prefactors = make_builder(&host_data.nuclear_form_prefactor); + prefactors.reserve(materials.num_isotopes()); + for (auto iso_id : range(IsotopeId{materials.num_isotopes()})) + { + prefactors.push_back( + this->calc_nuclear_form_prefactor(materials.get(iso_id))); + } + CELER_ENSURE(host_data.nuclear_form_prefactor.size() + == materials.num_isotopes()); + + // Build material data + if (host_data.params.is_combined) + { + std::vector inv_mass_cbrt_sq(materials.num_materials(), 0); + for (auto mat_id : range(MaterialId(materials.num_materials()))) + { + auto mat = materials.get(mat_id); + for (auto elcomp_id : range(ElementComponentId(mat.num_elements()))) + { + auto const& el_comp = mat.elements()[elcomp_id.get()]; + auto atomic_mass + = mat.make_element_view(elcomp_id).atomic_mass(); + inv_mass_cbrt_sq[mat_id.get()] + += el_comp.fraction + / std::pow(atomic_mass.value(), real_type(2) / 3); + } + inv_mass_cbrt_sq[mat_id.get()] *= mat.number_density(); + } + make_builder(&host_data.inv_mass_cbrt_sq) + .insert_back(inv_mass_cbrt_sq.begin(), inv_mass_cbrt_sq.end()); + } +} + +//---------------------------------------------------------------------------// +/*! + * Calculate the constant prefactors of the squared momentum transfer. + * + * This factor is used in the exponential and Gaussian nuclear form factor + * types: see Eqs. 2.262--2.264 of [LR15]. + * + * Specifically, it calculates \f$ (r_n/\bar h)^2 / 12 \f$, where \f$ r_n = + * 1.27 A^{0.27} \f$ [fm] is the nuclear radius. A special case is inherited + * from Geant for hydrogen targets. + */ +real_type WentzelOKVIParams::calc_nuclear_form_prefactor(IsotopeView const& iso) +{ + if (iso.atomic_number().get() == 1) + { + // TODO: Geant4 hardcodes a different prefactor for hydrogen + return real_type{1.5485e-6}; + } + + // The ratio has units of (MeV/c)^-2, so it's easier to convert the + // inverse which has units of MomentumSq, then invert afterwards + constexpr real_type ratio + = 1 + / native_value_to( + 12 + * ipow<2>(constants::hbar_planck + / (real_type(1.27) * units::femtometer))) + .value(); + return ratio + * fastpow(real_type(iso.atomic_mass_number().get()), + 2 * real_type(0.27)); +} + +//---------------------------------------------------------------------------// +/*! + * Get interpolated Mott coefficients for an element. + * + * These coefficients are used by the Lijian, Quing, Zhengming + * expression [PRM 8.48] for atomic numbers 1 <= Z <= 92. + * This data was taken from Geant4's G4MottData.hh file. + * + * For higher Z values, the PRM cites a numerical solution by Boschini + * et al (2013), but does not implement it. + */ +MottElementData::MottCoeffMatrix +WentzelOKVIParams::get_mott_coeff_matrix(AtomicNumber z) +{ + CELER_EXPECT(z); + + using CoeffMat = MottElementData::MottCoeffMatrix; + // clang-format off + static CoeffMat const mott_coeffs[] = { + //H....................................... + CoeffMat{{ + {1.e0,2.67363e-8,7.1153e-8,-9.7703e-8,-6.69132e-7,-3.09263e-7}, + {1.17182e-2,1.62222e-2,-5.90397e-5,-1.05585e-4,4.17873e-4,9.13843e-4}, + {-2.65955e-1,-7.29531e-1,-4.99796e-1,2.83507e-4,-9.09042e-4,-2.20244e-3}, + {-1.82348e-4,-8.86355e-5,-1.90554e-4,-2.49708e-4,6.35004e-4,1.73523e-3}, + {4.70966e-5,-4.09705e-6,3.75218e-5,8.05645e-5,-1.90534e-4,-5.42847e-4} + }}, + //He....................................... + CoeffMat{{ + {1.e0,3.76476e-8,-3.05313e-7,-3.27422e-7,2.44235e-6,4.08754e-6}, + {2.35767e-2,3.24642e-2,-6.37269e-4,-7.6916e-4,5.28004e-3,9.45642e-3}, + {-2.73743e-1,-7.40767e-1,-4.98195e-1,1.74337e-3,-1.25798e-2,-2.24046e-2}, + {-7.79128e-4,-4.14495e-4,-1.62657e-3,-1.37286e-3,1.04319e-2,1.83488e-2}, + {2.02855e-4,1.94598e-6,4.30102e-4,4.3218e-4,-3.31526e-3,-5.81788e-3} + }}, + //Li.................................................. + CoeffMat{{ + {1.e0,7.00357e-8,-3.15076e-7,-4.24915e-7,2.45516e-6,4.90187e-6}, + {3.55657e-2,4.87956e-2,-1.95525e-3,-2.7866e-3,1.6549e-2,3.11496e-2}, + {-2.81171e-1,-7.52015e-1,-4.95329e-1,5.83548e-3,-3.3983e-2,-6.55379e-2}, + {-1.83452e-3,-8.12746e-4,-3.84675e-3,-4.44467e-3,2.55871e-2,4.99483e-2}, + {4.79031e-4,-3.89615e-5,1.01022e-3,1.39133e-3,-7.99398e-3,-1.56366e-2} + }}, + //Be.................................................. + CoeffMat{{ + {1.e0,7.58881e-8,4.705e-8,2.48041e-7,-2.06053e-6,-1.97319e-6}, + {4.76788e-2,6.522e-2,-4.54331e-3,-6.50318e-3,3.76564e-2,7.17176e-2}, + {-2.88203e-1,-7.63217e-1,-4.90337e-1,1.22839e-2,-6.86398e-2,-1.35769e-1}, + {-3.37733e-3,-1.36514e-3,-7.51614e-3,-8.78592e-3,4.78572e-2,9.69021e-2}, + {8.81822e-4,-1.02577e-4,1.99797e-3,2.72661e-3,-1.48296e-2,-3.0106e-2} + }}, + //B.................................................... + CoeffMat{{ + {9.99999e-1,7.91498e-8,1.84164e-6,2.68534e-6,-1.8163e-5,-2.69021e-5}, + {5.98818e-2,8.17654e-2,-7.70811e-3,-1.12378e-2,6.38329e-2,1.25339e-1}, + {-2.94716e-1,-7.74405e-1,-4.8622e-1,1.77367e-2,-9.46825e-2,-2.01789e-1}, + {-5.52375e-3,-2.05348e-3,-9.44915e-3,-1.08135e-2,5.41024e-2,1.25257e-1}, + {1.44555e-3,-1.99404e-4,2.36742e-3,3.29655e-3,-1.64122e-2,-3.8375e-2} + }}, + //C................................................... + CoeffMat{{ + {9.99999e-1,7.68158e-8,5.18185e-6,7.34245e-6,-4.9478e-5,-7.71923e-5}, + {7.21461e-2,9.84618e-2,-1.06535e-2,-1.62358e-2,8.59238e-2,1.78727e-1}, + {-3.00622e-1,-7.85616e-1,-4.85735e-1,1.91563e-2,-8.10204e-2,-2.15386e-1}, + {-8.34809e-3,-2.85241e-3,-7.03252e-3,-7.56786e-3,1.44975e-2,8.79093e-2}, + {2.18964e-3,-3.42022e-4,1.3293e-3,2.20108e-3,-3.57927e-3,-2.5928e-2} + }}, + //N.................................................... + CoeffMat{{ + {9.99999e-1,8.36312e-8,1.09116e-5,1.47812e-5,-1.02733e-4,-1.62724e-4}, + {8.44142e-2,1.1531e-1,-1.1723e-2,-1.94732e-2,8.92604e-2,2.09303e-1}, + {-3.05743e-1,-7.96809e-1,-4.93957e-1,1.01607e-2,1.67761e-2,-1.05909e-1}, + {-1.2009e-2,-3.80678e-3,4.51195e-3,6.93472e-3,-1.12405e-1,-8.15484e-2}, + {3.16048e-3,-5.22237e-4,-2.58261e-3,-2.38303e-3,3.63393e-2,2.75127e-2} + }}, + //O................................................... + CoeffMat{{ + {9.99998e-1,1.57323e-8,1.77595e-5,2.56082e-5,-1.67537e-4,-2.73755e-4}, + {9.66438e-2,1.32264e-1,-9.53841e-3,-1.83707e-2,6.01664e-2,1.93357e-1}, + {-3.09969e-1,-8.0779e-1,-5.14392e-1,-1.67153e-2,2.3387e-1,1.916e-1}, + {-1.65906e-2,-5.11585e-3,2.80424e-2,3.94663e-2,-3.5572e-1,-4.39251e-1}, + {4.37866e-3,-6.81795e-4,-1.0152e-2,-1.24875e-2,1.11484e-1,1.38105e-1} + }}, + //F................................................. + CoeffMat{{ + {9.99997e-1,-8.06132e-8,2.49797e-5,3.8512e-5,-2.37451e-4,-3.99607e-4}, + {1.08782e-1,1.49306e-1,-2.50975e-3,-1.05471e-2,-1.64831e-2,1.05733e-1}, + {-3.13165e-1,-8.18489e-1,-5.50832e-1,-6.74447e-2,6.06357e-1,7.39717e-1}, + {-2.21976e-2,-6.84023e-3,6.65411e-2,9.48702e-2,-7.43989e-1,-1.03582e0}, + {5.87088e-3,-8.1048e-4,-2.21731e-2,-2.94422e-2,2.2954e-1,3.19669e-1} + }}, + //Ne................................................. + CoeffMat{{ + {9.99997e-1,-1.87404e-7,3.10276e-5,5.20007e-5,-2.98132e-4,-5.19259e-4}, + {1.20783e-1,1.66407e-1,1.06608e-2,6.48772e-3,-1.53031e-1,-7.59354e-2}, + {-3.15222e-1,-8.28793e-1,-6.0574e-1,-1.47812e-1,1.1576e0,1.58565e0}, + {-2.89055e-2,-9.08096e-3,1.21467e-1,1.77575e-1,-1.2911e0,-1.90333e0}, + {7.65342e-3,-8.85417e-4,-3.89092e-2,-5.4404e-2,3.93087e-1,5.79439e-1} + }}, + //Na................................................. + CoeffMat{{ + {9.99996e-1,-2.44548e-7,3.31019e-5,6.29483e-5,-3.24667e-4,-5.95527e-4}, + {1.32615e-1,1.83566e-1,3.04158e-2,3.40925e-2,-3.54681e-1,-3.63044e-1}, + {-3.16092e-1,-8.38704e-1,-6.78558e-1,-2.59346e-1,1.88547e0,2.73632e0}, + {-3.67233e-2,-1.18139e-2,1.91089e-1,2.87408e-1,-1.98397e0,-3.03075e0}, + {9.72033e-3,-9.2638e-4,-5.95654e-2,-8.69829e-2,5.95744e-1,9.10242e-1} + }}, + //Mg................................................. + CoeffMat{{ + {9.99995e-1,-2.12227e-7,2.95645e-5,6.92848e-5,-3.02153e-4,-6.05145e-4}, + {1.44258e-1,2.00775e-1,5.67845e-2,7.35166e-2,-6.22861e-1,-7.62213e-1}, + {-3.15754e-1,-8.48196e-1,-7.67318e-1,-4.02984e-1,2.77477e0,4.18114e0}, + {-4.56307e-2,-1.50425e-2,2.72232e-1,4.23528e-1,-2.79606e0,-4.38863e0}, + {1.2056e-2,-9.44637e-4,-8.28738e-2,-1.26564e-1,8.26726e-1,1.29882e0} + }}, + //Al..................................................... + CoeffMat{{ + {9.99995e-1,-4.03407e-8,1.86047e-5,6.85201e-5,-2.14503e-4,-5.22528e-4}, + {1.55704e-1,2.18048e-1,8.88994e-2,1.24878e-1,-9.51331e-1,-1.26824e0}, + {-3.14244e-1,-8.57322e-1,-8.6719e-1,-5.75787e-1,3.78571e0,5.87052e0}, + {-5.55526e-2,-1.86861e-2,3.5886e-1,5.81094e-1,-3.67623e0,-5.908e0}, + {1.46269e-2,-9.79742e-4,-1.06652e-1,-1.71226e-1,1.06737e0,1.71918e0} + }}, + //Si..................................................... + CoeffMat{{ + {9.99994e-1,3.00267e-7,-1.1184e-6,5.88256e-5,-4.78456e-5,-3.25731e-4}, + {1.6696e-1,2.35405e-1,1.25215e-1,1.87646e-1,-1.32685e0,-1.86549e0}, + {-3.1163e-1,-8.66152e-1,-9.71254e-1,-7.72715e-1,4.85654e0,7.7215e0}, + {-6.63778e-2,-2.26481e-2,4.42898e-1,7.53182e-1,-4.55172e0,-7.4867e0}, + {1.73883e-2,-1.07669e-3,-1.28075e-1,-2.18389e-1,1.29217e0,2.13475e0} + }}, + //P..................................................... + CoeffMat{{ + {9.99994e-1,8.94829e-7,-2.98434e-5,3.82193e-5,2.00584e-4,-6.40482e-6}, + {1.78039e-1,2.52912e-1,1.63761e-1,2.60132e-1,-1.73287e0,-2.53185e0}, + {-3.08007e-1,-8.74905e-1,-1.07146e0,-9.85062e-1,5.91697e0,9.63265e0}, + {-7.79747e-2,-2.66797e-2,5.15288e-1,9.29261e-1,-5.34252e0,-9.00574e0}, + {2.02892e-2,-1.33011e-3,-1.44039e-1,-2.6433e-1,1.4736e0,2.50398e0} + }}, + //S.................................................... + CoeffMat{{ + {9.99994e-1,1.75397e-6,-6.7331e-5,6.29524e-6,5.29623e-4,4.35288e-4}, + {1.88968e-1,2.70612e-1,2.01975e-1,3.40574e-1,-2.14737e0,-3.23836e0}, + {-3.03499e-1,-8.83717e-1,-1.15816e0,-1.20414e0,6.88176e0,1.14841e1}, + {-9.01806e-2,-3.06202e-2,5.65581e-1,1.09902e0,-5.95552e0,-1.03302e1}, + {2.32694e-2,-1.80614e-3,-1.51041e-1,-3.05449e-1,1.58037e0,2.78083e0} + }}, + //Cl.................................................... + CoeffMat{{ + {9.99994e-1,3.07931e-6,-1.11876e-4,-4.10164e-5,9.17095e-4,9.80145e-4}, + {1.99765e-1,2.88611e-1,2.37501e-1,4.25803e-1,-2.55105e0,-3.95585e0}, + {-2.98206e-1,-8.9293e-1,-1.22279e0,-1.4169e0,7.67836e0,1.31601e1}, + {-1.02865e-1,-3.40967e-2,5.84677e-1,1.24786e0,-6.31301e0,-1.13328e1}, + {2.628e-2,-2.63995e-3,-1.46076e-1,-3.36795e-1,1.58677e0,2.92251e0} + }}, + //Ar.................................................... + CoeffMat{{ + {9.99993e-1,4.49776e-6,-1.65136e-4,-9.76754e-5,1.39664e-3,1.66293e-3}, + {2.10469e-1,3.06924e-1,2.66793e-1,5.13797e-1,-2.90958e0,-4.63816e0}, + {-2.92294e-1,-9.0256e-1,-1.25307e0,-1.6147e0,8.18574e0,1.44912e1}, + {-1.15831e-1,-3.70891e-2,5.59807e-1,1.36619e0,-6.28824e0,-1.18327e1}, + {2.92513e-2,-3.84903e-3,-1.24976e-1,-3.55149e-1,1.45127e0,2.86925e0} + }}, + //K................................................. + CoeffMat{{ + {9.99993e-1,6.01488e-6,-2.22125e-4,-1.61774e-4,1.92058e-3,2.41975e-3}, + {2.21091e-1,3.25648e-1,2.87732e-1,6.01632e-1,-3.20436e0,-5.25724e0}, + {-2.85814e-1,-9.12889e-1,-1.24213e0,-1.78635e0,8.34196e0,1.53776e1}, + {-1.29005e-1,-3.92986e-2,4.84255e-1,1.44206e0,-5.81999e0,-1.17275e1}, + {3.21555e-2,-5.54272e-3,-8.56572e-2,-3.56589e-1,1.1548e0,2.58829e0} + }}, + //Ca................................................. + CoeffMat{{ + {9.99993e-1,8.01467e-6,-2.79242e-4,-2.3682e-4,2.45459e-3,3.21683e-3}, + {2.31651e-1,3.44948e-1,2.9782e-1,6.85187e-1,-3.41294e0,-5.77715e0}, + {-2.78858e-1,-9.24428e-1,-1.18215e0,-1.91691e0,8.07489e0,1.56969e1}, + {-1.42276e-1,-4.01888e-2,3.50466e-1,1.45983e0,-4.83806e0,-1.08936e1}, + {3.49529e-2,-7.90933e-3,-2.58002e-2,-3.36028e-1,6.7574e-1,2.04052e0} + }}, + //Sc................................................. + CoeffMat{{ + {9.99992e-1,1.04277e-5,-3.35126e-4,-3.21042e-4,2.98507e-3,4.03325e-3}, + {2.42172e-1,3.64954e-1,2.94606e-1,7.60693e-1,-3.51409e0,-6.1646e0}, + {-2.71512e-1,-9.37543e-1,-1.0657e0,-1.99328e0,7.31863e0,1.53396e1}, + {-1.5554e-1,-3.93862e-2,1.51477e-1,1.40614e0,-3.28024e0,-9.22338e0}, + {3.76066e-2,-1.10812e-2,5.66831e-2,-2.8921e-1,-4.60274e-3,1.19273e0} + }}, + //Ti................................................. + CoeffMat{{ + {9.99992e-1,1.30838e-5,-3.8407e-4,-4.09294e-4,3.47025e-3,4.81071e-3}, + {2.52646e-1,3.85718e-1,2.7633e-1,8.25665e-1,-3.48888e0,-6.39017e0}, + {-2.63758e-1,-9.52326e-1,-8.88179e-1,-2.0072e0,6.0196e0,1.42157e1}, + {-1.68806e-1,-3.68095e-2,-1.16429e-1,1.27309e0,-1.09991e0,-6.63409e0}, + {4.01184e-2,-1.50938e-2,1.6274e-1,-2.1376e-1,-8.99142e-1,2.08129e-2} + }}, + //V................................................. + CoeffMat{{ + {9.99991e-1,1.59363e-5,-4.27315e-4,-5.01461e-4,3.91365e-3,5.55096e-3}, + {2.63096e-1,4.07357e-1,2.40649e-1,8.7645e-1,-3.31641e0,-6.42069e0}, + {-2.55682e-1,-9.6909e-1,-6.43149e-1,-1.94678e0,4.11915e0,1.22255e1}, + {-1.81974e-1,-3.21462e-2,-4.58817e-1,1.04913e0,1.75388e0,-3.03434e0}, + {4.24541e-2,-2.00595e-2,2.93916e-1,-1.06138e-1,-2.02203e0,-1.50181e0} + }}, + //Cr................................................. + CoeffMat{{ + {9.9999e-1,1.8895e-5,-4.59994e-4,-5.93663e-4,4.27684e-3,6.2009e-3}, + {2.73504e-1,4.2999e-1,1.86473e-1,9.09921e-1,-2.98441e0,-6.23344e0}, + {-2.47225e-1,-9.88118e-1,-3.28759e-1,-1.80252e0,1.5909e0,9.30968e0}, + {-1.951e-1,-2.51242e-2,-8.76244e-1,7.25592e-1,5.2965e0,1.62177e0}, + {4.46307e-2,-2.60754e-2,4.50054e-1,3.61638e-2,-3.37524e0,-3.38601e0} + }}, + //Mn................................................. + CoeffMat{{ + {9.99989e-1,2.18906e-5,-4.79199e-4,-6.83164e-4,4.53467e-3,6.72485e-3}, + {2.83854e-1,4.53722e-1,1.12877e-1,9.23203e-1,-2.48247e0,-5.80855e0}, + {-2.38338e-1,-1.00965e0,5.61453e-2,-1.56605e0,-1.58353e0,5.42138e0}, + {-2.08226e-1,-1.55138e-2,-1.36846e0,2.95168e-1,9.53392e0,7.3655e0}, + {4.66619e-2,-3.32255e-2,6.30711e-1,2.15167e-1,-4.95748e0,-5.63744e0} + }}, + //Fe................................................... + CoeffMat{{ + {9.99987e-1,2.482e-5,-4.82895e-4,-7.67488e-4,4.669e-3,7.09581e-3}, + {2.94123e-1,4.78653e-1,1.93256e-2,9.13569e-1,-1.80354e0,-5.1304e0}, + {-2.28945e-1,-1.0339e0,5.1121e-1,-1.22996e0,-5.40939e0,5.3149e-1}, + {-2.21426e-1,-3.12679e-3,-1.93353e0,-2.48179e-1,1.44576e1,1.42077e1}, + {4.85718e-2,-4.15796e-2,8.34879e-1,4.32418e-1,-6.76244e0,-8.25464e0} + }}, + //Co.................................................. + CoeffMat{{ + {9.99985e-1,2.76168e-5,-4.67522e-4,-8.43978e-4,4.65008e-3,7.27476e-3}, + {3.0428e-1,5.04856e-1,-9.42256e-2,8.78905e-1,-9.44256e-1,-4.18847e0}, + {-2.18939e-1,-1.06098e0,1.03431e0,-7.89123e-1,-9.87815e0,-5.36997e0}, + {-2.34805e-1,1.21281e-2,-2.56765e0,-9.08015e-1,2.00436e1,2.21379e1}, + {5.03942e-2,-5.11768e-2,1.0609e0,6.88631e-1,-8.77867e0,-1.12289e1} + }}, + //Ni.................................................. + CoeffMat{{ + {9.99982e-1,3.00792e-5,-4.33447e-4,-9.09366e-4,4.47786e-3,7.25072e-3}, + {3.14283e-1,5.32444e-1,-2.27528e-1,8.16951e-1,9.53704e-2,-2.9773e0}, + {-2.08186e-1,-1.09112e0,1.62237e0,-2.38394e-1,-1.49699e1,-1.22743e1}, + {-2.48493e-1,3.04543e-2,-3.26601e0,-1.6876e0,2.62568e1,3.11259e1}, + {5.21709e-2,-6.20908e-2,1.30683e0,9.84353e-1,-1.09909e1,-1.45449e1} + }}, + //Cu................................................... + CoeffMat{{ + {9.99979e-1,3.24569e-5,-3.76717e-4,-9.64e-4,4.11997e-3,6.98761e-3}, + {3.24082e-1,5.61534e-1,-3.79993e-1,7.25313e-1,1.31342e0,-1.49163e0}, + {-1.96521e-1,-1.12457e0,2.27091e0,4.27678e-1,-2.06558e1,-2.0169e1}, + {-2.62655e-1,5.20812e-2,-4.02224e0,-2.59048e0,3.30508e1,4.1135e1}, + {5.39556e-2,-7.44077e-2,1.57015e0,1.32022e0,-1.33802e1,-1.81848e1} + }}, + //Zn................................................... + CoeffMat{{ + {9.99976e-1,3.39628e-5,-2.98845e-4,-9.99651e-4,3.58523e-3,6.47782e-3}, + {3.33647e-1,5.92068e-1,-5.5102e-1,6.0363e-1,2.70712e0,2.67853e-1}, + {-1.83849e-1,-1.16095e0,2.97559e0,1.20728e0,-2.69053e1,-2.90214e1}, + {-2.77375e-1,7.65838e-2,-4.83023e0,-3.61238e0,4.03786e1,5.21084e1}, + {5.57745e-2,-8.79952e-2,1.84848e0,1.69425e0,-1.59274e1,-2.21242e1} + }}, + //Ga................................................... + CoeffMat{{ + {9.99972e-1,3.48473e-5,-1.97828e-4,-1.01659e-3,2.85509e-3,5.69725e-3}, + {3.42904e-1,6.24164e-1,-7.39023e-1,4.50306e-1,4.26553e0,2.29299e0}, + {-1.6994e-1,-1.20051e0,3.72868e0,2.10276e0,-3.36594e1,-3.87726e1}, + {-2.92882e-1,1.04194e-1,-5.68033e0,-4.75339e0,4.81633e1,6.39614e1}, + {5.77008e-2,-1.02941e-1,2.13826e0,2.10589e0,-1.86034e1,-2.63294e1} + }}, + //Ge................................................. + CoeffMat{{ + {9.99968e-1,3.47804e-5,-7.11898e-5,-1.01028e-3,1.90919e-3,4.61426e-3}, + {3.51793e-1,6.57815e-1,-9.42175e-1,2.6494e-1,5.97606e0,4.57222e0}, + {-1.54595e-1,-1.24305e0,4.5217e0,3.11205e0,-4.08539e1,-4.93518e1}, + {-3.09367e-1,1.34661e-1,-6.56214e0,-6.00882e0,5.63228e1,7.65965e1}, + {5.97948e-2,-1.19174e-1,2.4357e0,2.553e0,-2.13779e1,-3.07628e1} + }}, + //As................................................. + CoeffMat{{ + {9.99963e-1,3.37519e-5,8.13037e-5,-9.78638e-4,7.41412e-4,3.21498e-3}, + {3.60246e-1,6.93065e-1,-1.1588e0,4.68519e-2,7.82662e0,7.09456e0}, + {-1.37615e-1,-1.28855e0,5.34673e0,4.23402e0,-4.84263e1,-6.06897e1}, + {-3.27017e-1,1.67941e-1,-7.46593e0,-7.37492e0,6.47773e1,8.99181e1}, + {6.21158e-2,-1.36692e-1,2.73726e0,3.03374e0,-2.42209e1,-3.53873e1} + }}, + //Se................................................. + CoeffMat{{ + {9.99958e-1,3.14983e-5,2.59741e-4,-9.19008e-4,-6.49202e-4,1.49008e-3}, + {3.68196e-1,7.29888e-1,-1.38664e0,-2.03807e-1,9.80062e0,9.84155e0}, + {-1.18788e-1,-1.33674e0,6.194e0,5.46446e0,-5.62994e1,-7.26915e1}, + {-3.46033e-1,2.03729e-1,-8.38012e0,-8.84466e0,7.34319e1,1.03805e2}, + {6.47255e-2,-1.55409e-1,3.03879e0,3.54516e0,-2.7098e1,-4.01572e1} + }}, + //Br................................................. + CoeffMat{{ + {9.99952e-1,2.79961e-5,4.60479e-4,-8.33486e-4,-2.23214e-3,-5.16285e-4}, + {3.7558e-1,7.68292e-1,-1.62392e0,-4.87674e-1,1.18854e1,1.28025e1}, + {-9.79355e-2,-1.38749e0,7.0556e0,6.80198e0,-6.44113e1,-8.52915e1}, + {-3.66572e-1,2.41863e-1,-9.29524e0,-1.04141e1,8.22093e1,1.18166e2}, + {6.76714e-2,-1.75289e-1,3.33691e0,4.08538e0,-2.99809e1,-4.50381e1} + }}, + //Kr................................................. + CoeffMat{{ + {9.99947e-1,2.2952e-5,6.82639e-4,-7.17139e-4,-4.00522e-3,-2.81601e-3}, + {3.82332e-1,8.08194e-1,-1.86848e0,-8.03415e-1,1.4064e1,1.5956e1}, + {-7.48735e-2,-1.44031e0,7.92233e0,8.2382e0,-7.26858e1,-9.83865e1}, + {-3.88797e-1,2.81837e-1,-1.02005e1,-1.20716e1,9.10175e1,1.32873e2}, + {7.10017e-2,-1.96182e-1,3.6278e0,4.65e0,-3.28365e1,-4.99823e1} + }}, + //Rb................................................. + CoeffMat{{ + {9.99941e-1,1.63607e-5,9.28242e-4,-5.69257e-4,-5.98245e-3,-5.42476e-3}, + {3.88363e-1,8.49548e-1,-2.11701e0,-1.15003e0,1.63119e1,1.92737e1}, + {-4.93364e-2,-1.49491e0,8.78129e0,9.76596e0,-8.10212e1,-1.11851e2}, + {-4.12953e-1,3.23322e-1,-1.10814e1,-1.38073e1,9.97386e1,1.47775e2}, + {7.47911e-2,-2.18001e-1,3.90642e0,5.23509e0,-3.56231e1,-5.49352e1} + }}, + //Sr................................................. + CoeffMat{{ + {9.99935e-1,8.3152e-6,1.19244e-3,-3.95258e-4,-8.12525e-3,-8.28836e-3}, + {3.93603e-1,8.92358e-1,-2.36688e0,-1.52786e0,1.86071e1,2.27314e1}, + {-2.11345e-2,-1.55112e0,9.62205e0,1.13827e1,-8.93266e1,-1.25574e2}, + {-4.39199e-1,3.66161e-1,-1.19262e1,-1.56159e1,1.08267e2,1.62736e2}, + {7.90853e-2,-2.40719e-1,4.16872e0,5.83839e0,-3.8304e1,-5.98484e1} + }}, + //Y................................................. + CoeffMat{{ + {9.99929e-1,-1.65608e-6,1.47476e-3,-1.8625e-4,-1.04363e-2,-1.14318e-2}, + {3.97989e-1,9.3643e-1,-2.61569e0,-1.93387e0,2.09305e1,2.63021e1}, + {9.88744e-3,-1.60813e0,1.04351e1,1.3074e1,-9.75207e1,-1.39437e2}, + {-4.67657e-1,4.09494e-1,-1.2724e1,-1.74797e1,1.16508e2,1.77614e2}, + {8.39169e-2,-2.64077e-1,4.41092e0,6.45344e0,-4.08455e1,-6.46698e1} + }}, + //Zr................................................. + CoeffMat{{ + {9.99922e-1,-1.37624e-5,1.77335e-3,5.994e-5,-1.29013e-2,-1.48443e-2}, + {4.0145e-1,9.816e-1,-2.86053e0,-2.36562e0,2.32591e1,2.99551e1}, + {4.3915e-2,-1.66522e0,1.12093e1,1.48279e1,-1.05512e2,-1.5331e2}, + {-4.98475e-1,4.52607e-1,-1.34628e1,-1.93836e1,1.24357e2,1.92256e2}, + {8.93258e-2,-2.87869e-1,4.62891e0,7.07473e0,-4.32116e1,-6.93461e1} + }}, + //Nb................................................. + CoeffMat{{ + {9.99916e-1,-2.78975e-5,2.08753e-3,3.41607e-4,-1.55141e-2,-1.85156e-2}, + {4.03889e-1,1.02783e0,-3.09762e0,-2.82183e0,2.55629e1,3.36548e1}, + {8.12122e-2,-1.7221e0,1.1931e1,1.66362e1,-1.13185e2,-1.67048e2}, + {-5.3188e-1,4.95236e-1,-1.41275e1,-2.13166e1,1.31687e2,2.06495e2}, + {9.53781e-2,-3.12041e-1,4.81765e0,7.69813e0,-4.53587e1,-7.38182e1} + }}, + //Mo................................................. + CoeffMat{{ + {9.9991e-1,-4.3988e-5,2.41439e-3,6.54636e-4,-1.82489e-2,-2.24062e-2}, + {4.05249e-1,1.07495e0,-3.32445e0,-3.30073e0,2.78221e1,3.7376e1}, + {1.21896e-1,-1.77813e0,1.25907e1,1.84893e1,-1.20461e2,-1.80542e2}, + {-5.67942e-1,5.36741e-1,-1.4708e1,-2.32664e1,1.38408e2,2.20203e2}, + {1.02086e-1,-3.36419e-1,4.97373e0,8.31905e0,-4.72562e1,-7.80409e1} + }}, + //Tc................................................. + CoeffMat{{ + {9.99904e-1,-6.22073e-5,2.75066e-3,1.00063e-3,-2.10817e-2,-2.64933e-2}, + {4.05472e-1,1.12279e0,-3.53825e0,-3.79964e0,3.00132e1,4.10871e1}, + {1.66089e-1,-1.83255e0,1.31785e1,2.03743e1,-1.27252e2,-1.93664e2}, + {-6.06729e-1,5.76392e-1,-1.51936e1,-2.52171e1,1.44423e2,2.33232e2}, + {1.09463e-1,-3.60803e-1,5.09358e0,8.93178e0,-4.88706e1,-8.19632e1} + }}, + //Ru................................................. + CoeffMat{{ + {9.99898e-1,-8.26232e-5,3.09353e-3,1.37924e-3,-2.39902e-2,-3.07502e-2}, + {4.04498e-1,1.17115e0,-3.73628e0,-4.3158e0,3.21142e1,4.4758e1}, + {2.13904e-1,-1.88458e0,1.36848e1,2.22782e1,-1.33471e2,-2.06293e2}, + {-6.48301e-1,6.13439e-1,-1.55742e1,-2.71531e1,1.49638e2,2.45444e2}, + {1.17516e-1,-3.8499e-1,5.17383e0,9.5307e0,-5.01709e1,-8.55368e1} + }}, + //Rh................................................. + CoeffMat{{ + {9.99893e-1,-1.05293e-4,3.43932e-3,1.78981e-3,-2.69462e-2,-3.51436e-2}, + {4.02272e-1,1.21982e0,-3.91575e0,-4.84621e0,3.41013e1,4.83563e1}, + {2.65438e-1,-1.93341e0,1.40998e1,2.41876e1,-1.39034e2,-2.183e2}, + {-6.92691e-1,6.47135e-1,-1.58398e1,-2.90581e1,1.5396e2,2.56697e2}, + {1.26243e-1,-4.08783e-1,5.21125e0,1.011e1,-5.11262e1,-8.8713e1} + }}, + //Pd................................................. + CoeffMat{{ + {9.99888e-1,-1.30188e-4,3.78436e-3,2.22971e-3,-2.99185e-2,-3.96315e-2}, + {3.9875e-1,1.26855e0,-4.07427e0,-5.38796e0,3.5955e1,5.18546e1}, + {3.20748e-1,-1.97817e0,1.44153e1,2.6089e1,-1.43866e2,-2.29578e2}, + {-7.39891e-1,6.76679e-1,-1.59819e1,-3.09162e1,1.57312e2,2.66867e2}, + {1.3563e-1,-4.31972e-1,5.2031e0,1.06642e1,-5.17104e1,-9.14493e1} + }}, + //Ag................................................. + CoeffMat{{ + {9.99884e-1,-1.57459e-4,4.12052e-3,2.69984e-3,-3.28468e-2,-4.41512e-2}, + {3.93891e-1,1.3171e0,-4.20944e0,-5.93691e0,3.76535e1,5.52188e1}, + {3.79857e-1,-2.01791e0,1.46238e1,2.79653e1,-1.47893e2,-2.39997e2}, + {-7.89852e-1,7.01196e-1,-1.59929e1,-3.27077e1,1.59613e2,2.75813e2}, + {1.45646e-1,-4.54328e-1,5.14701e0,1.11863e1,-5.18977e1,-9.36987e1} + }}, + //Cd................................................. + CoeffMat{{ + {9.99881e-1,-1.86816e-4,4.45491e-3,3.19574e-3,-3.57812e-2,-4.87457e-2}, + {3.87612e-1,1.36524e0,-4.31686e0,-6.49005e0,3.91612e1,5.84032e1}, + {4.4294e-1,-2.05189e0,1.47103e1,2.98033e1,-1.50988e2,-2.49389e2}, + {-8.42683e-1,7.20045e-1,-1.58579e1,-3.44169e1,1.60733e2,2.83354e2}, + {1.56312e-1,-4.75705e-1,5.0381e0,1.16708e1,-5.16454e1,-9.53999e1} + }}, + //In................................................. + CoeffMat{{ + {9.99878e-1,-2.1848e-4,4.78376e-3,3.71897e-3,-3.86926e-2,-5.33866e-2}, + {3.79879e-1,1.41265e0,-4.39391e0,-7.04282e0,4.0455e1,6.13721e1}, + {5.10001e-1,-2.07896e0,1.46666e1,3.15843e1,-1.53072e2,-2.57624e2}, + {-8.98306e-1,7.32213e-1,-1.55686e1,-3.60229e1,1.60591e2,2.89348e2}, + {1.67591e-1,-4.95837e-1,4.87389e0,1.21106e1,-5.09274e1,-9.6506e1} + }}, + //Sn................................................ + CoeffMat{{ + {9.99876e-1,-2.52173e-4,5.09558e-3,4.26518e-3,-4.14939e-2,-5.79712e-2}, + {3.70689e-1,1.45908e0,-4.43983e0,-7.59192e0,4.15261e1,6.4108e1}, + {5.80924e-1,-2.09831e0,1.44911e1,3.32945e1,-1.54118e2,-2.64634e2}, + {-9.56518e-1,7.37009e-1,-1.51243e1,-3.75098e1,1.59159e2,2.93723e2}, + {1.79397e-1,-5.14573e-1,4.65426e0,1.25003e1,-4.97362e1,-9.69932e1} + }}, + //Sb................................................. + CoeffMat{{ + {9.99874e-1,-2.87917e-4,5.39066e-3,4.83313e-3,-4.41852e-2,-6.24953e-2}, + {3.60001e-1,1.5042e0,-4.4515e0,-8.1327e0,4.23476e1,6.65709e1}, + {6.55716e-1,-2.10882e0,1.41741e1,3.49157e1,-1.54034e2,-2.70276e2}, + {-1.01724e0,7.33501e-1,-1.45156e1,-3.88572e1,1.56348e2,2.96328e2}, + {1.9169e-1,-5.3169e-1,4.37639e0,1.28328e1,-4.80435e1,-9.68125e1} + }}, + //Te................................................. + CoeffMat{{ + {9.99874e-1,-3.25847e-4,5.67083e-3,5.42945e-3,-4.67877e-2,-6.701e-2}, + {3.47783e-1,1.54767e0,-4.42606e0,-8.65995e0,4.28934e1,6.87183e1}, + {7.34353e-1,-2.10941e0,1.3707e1,3.64273e1,-1.52736e2,-2.74401e2}, + {-1.08036e0,7.20771e-1,-1.37344e1,-4.00427e1,1.52073e2,2.97007e2}, + {2.04417e-1,-5.46972e-1,4.03784e0,1.31007e1,-4.58224e1,-9.59123e1} + }}, + //I................................................. + CoeffMat{{ + {9.99875e-1,-3.65406e-4,5.92115e-3,6.03332e-3,-4.91674e-2,-7.12951e-2}, + {3.34062e-1,1.58924e0,-4.36361e0,-9.17241e0,4.31638e1,7.05522e1}, + {8.16604e-1,-2.09933e0,1.30915e1,3.78235e1,-1.50231e2,-2.77015e2}, + {-1.14554e0,6.98298e-1,-1.2784e1,-4.10591e1,1.46348e2,2.9577e2}, + {2.17451e-1,-5.60344e-1,3.63994e0,1.33015e1,-4.30802e1,-9.42981e1} + }}, + //Xe.................................................. + CoeffMat{{ + {9.99877e-1,-4.06637e-4,6.14278e-3,6.64937e-3,-5.13391e-2,-7.53881e-2}, + {3.18822e-1,1.62854e0,-4.26169e0,-9.66459e0,4.31357e1,7.20331e1}, + {9.02373e-1,-2.07743e0,1.23209e1,3.90832e1,-1.46445e2,-2.77984e2}, + {-1.21259e0,6.65174e-1,-1.16583e1,-4.18839e1,1.39104e2,2.92479e2}, + {2.3071e-1,-5.71608e-1,3.18104e0,1.34277e1,-3.97963e1,-9.19255e1} + }}, + //Cs.................................................. + CoeffMat{{ + {9.99881e-1,-4.49175e-4,6.32631e-3,7.27184e-3,-5.32297e-2,-7.91972e-2}, + {3.02086e-1,1.66529e0,-4.11996e0,-1.0133e1,4.28031e1,7.31464e1}, + {9.91429e-1,-2.04292e0,1.13959e1,4.01931e1,-1.41369e2,-2.7726e2}, + {-1.28117e0,6.2087e-1,-1.03596e1,-4.25024e1,1.30338e2,2.8709e2}, + {2.44065e-1,-5.80702e-1,2.66217e0,1.34745e1,-3.59721e1,-8.87824e1} + }}, + //Ba.................................................. + CoeffMat{{ + {9.99886e-1,-4.93136e-4,6.47914e-3,7.90029e-3,-5.48909e-2,-8.27736e-2}, + {2.8384e-1,1.69905e0,-3.93534e0,-1.05724e1,4.21392e1,7.38521e1}, + {1.08366e0,-1.99451e0,1.03076e1,4.1134e1,-1.34917e2,-2.7471e2}, + {-1.35107e0,5.64397e-1,-8.87993e0,-4.28945e1,1.19972e2,2.7947e2}, + {2.57431e-1,-5.87417e-1,2.08111e0,1.34353e1,-3.15844e1,-8.48268e1} + }}, + //La.................................................. + CoeffMat{{ + {9.99892e-1,-5.38122e-4,6.59206e-3,8.52848e-3,-5.62508e-2,-8.60251e-2}, + {2.6413e-1,1.72952e0,-3.70812e0,-1.09794e1,4.1143e1,7.4141e1}, + {1.17876e0,-1.93141e0,9.05928e0,4.18928e1,-1.27099e2,-2.70309e2}, + {-1.42185e0,4.95292e-1,-7.22427e0,-4.30461e1,1.08024e2,2.69601e2}, + {2.70647e-1,-5.91727e-1,1.43982e0,1.33056e1,-2.66423e1,-8.00557e1} + }}, + //Ce.................................................. + CoeffMat{{ + {9.999e-1,-5.83896e-4,6.66373e-3,9.15389e-3,-5.72988e-2,-8.89368e-2}, + {2.42961e-1,1.75635e0,-3.43619e0,-1.13482e1,3.97928e1,7.39747e1}, + {1.27654e0,-1.85263e0,7.64589e0,4.24492e1,-1.17853e2,-2.63933e2}, + {-1.49322e0,4.12926e-1,-5.38909e0,-4.29363e1,9.44419e1,2.57363e2}, + {2.83597e-1,-5.93557e-1,7.3754e-1,1.30785e1,-2.11316e1,-7.44329e1} + }}, + //Pr.................................................. + CoeffMat{{ + {9.99909e-1,-6.301e-4,6.69295e-3,9.77077e-3,-5.8022e-2,-9.14799e-2}, + {2.20373e-1,1.77919e0,-3.11907e0,-1.16754e1,3.80818e1,7.33384e1}, + {1.37668e0,-1.75734e0,6.06819e0,4.27905e1,-1.07169e2,-2.55543e2}, + {-1.56475e0,3.16885e-1,-3.37687e0,-4.25518e1,7.92269e1,2.42726e2}, + {2.96124e-1,-5.92915e-1,-2.45072e-2,1.27501e1,-1.50562e1,-6.79514e1} + }}, + //Nd.................................................. + CoeffMat{{ + {9.9992e-1,-6.77211e-4,6.67917e-3,1.03866e-2,-5.84208e-2,-9.36775e-2}, + {1.96406e-1,1.79753e0,-2.7553e0,-1.19533e1,3.59929e1,7.21943e1}, + {1.47887e0,-1.64416e0,4.32381e0,4.28906e1,-9.50034e1,-2.45021e2}, + {-1.636e0,2.06169e-1,-1.18704e0,-4.18661e1,6.23476e1,2.25579e2}, + {3.08067e-1,-5.89622e-1,-8.4608e-1,1.23116e1,-8.40939e0,-6.05788e1} + }}, + //Pm.................................................. + CoeffMat{{ + {9.99932e-1,-7.23543e-4,6.61026e-3,1.09756e-2,-5.83892e-2,-9.53462e-2}, + {1.71135e-1,1.81127e0,-2.34665e0,-1.21818e1,3.35368e1,7.05547e1}, + {1.58269e0,-1.5131e0,2.4214e0,4.27487e1,-8.14091e1,-2.32422e2}, + {-1.70642e0,8.13302e-2,1.16945e0,-4.0879e1,4.38728e1,2.05995e2}, + {3.19226e-1,-5.84023e-1,-1.72313e0,1.17636e1,-1.21776e0,-5.23434e1} + }}, + //Sm.................................................. + CoeffMat{{ + {9.99946e-1,-7.69368e-4,6.48923e-3,1.15425e-2,-5.79523e-2,-9.65284e-2}, + {1.44613e-1,1.81998e0,-1.89196e0,-1.23548e1,3.06998e1,6.83894e1}, + {1.68776e0,-1.36316e0,3.59801e-1,4.23443e1,-6.63559e1,-2.17658e2}, + {-1.77552e0,-5.81387e-2,3.69184e0,-3.95701e1,2.37855e1,1.83897e2}, + {3.2942e-1,-5.76119e-1,-2.65496e0,1.10995e1,6.52027e0,-4.32241e1} + }}, + //Eu.................................................. + CoeffMat{{ + {9.99962e-1,-8.13654e-4,6.31217e-3,1.20763e-2,-5.70782e-2,-9.71683e-2}, + {1.16915e-1,1.82343e0,-1.39208e0,-1.24689e1,2.74833e1,6.56914e1}, + {1.79364e0,-1.194e0,-1.85513e0,4.16665e1,-4.9868e1,-2.00724e2}, + {-1.84272e0,-2.12012e-1,6.37189e0,-3.79293e1,2.12742e0,1.59295e2}, + {3.38449e-1,-5.66167e-1,-3.63838e0,1.03165e1,1.47866e1,-3.32296e1} + }}, + //Gd.................................................. + CoeffMat{{ + {9.9998e-1,-8.57331e-4,6.07736e-3,1.25819e-2,-5.57471e-2,-9.72475e-2}, + {8.8154e-2,1.82105e0,-8.47718e-1,-1.25186e1,2.38891e1,6.2451e1}, + {1.89976e0,-1.00428e0,-4.21799e0,4.0697e1,-3.19703e1,-1.81604e2}, + {-1.90731e0,-3.8108e-1,9.2019e0,-3.59386e1,-2.10594e1,1.32192e2}, + {3.46059e-1,-5.54093e-1,-4.67039e0,9.40924e0,2.35631e1,-2.23653e1} + }}, + //Tb.................................................. + CoeffMat{{ + {9.99999e-1,-8.99343e-4,5.7868e-3,1.30512e-2,-5.39764e-2,-9.67723e-2}, + {5.84024e-2,1.81266e0,-2.59193e-1,-1.25002e1,1.99139e1,5.86545e1}, + {2.00567e0,-7.93885e-1,-6.72461e0,3.94244e1,-1.26714e1,-1.60273e2}, + {-1.96873e0,-5.64782e-1,1.21753e1,-3.3588e1,-4.57487e1,1.02578e2}, + {3.5205e-1,-5.40282e-1,-5.74837e0,8.37487e0,3.28364e1,-1.06341e1} + }}, + //Dy.................................................. + CoeffMat{{ + {1.00002e0,-9.39311e-4,5.43607e-3,1.34763e-2,-5.17249e-2,-9.56765e-2}, + {2.78005e-2,1.79795e0,3.70737e-1,-1.24114e1,1.55759e1,5.43184e1}, + {2.1107e0,-5.62473e-1,-9.36239e0,3.78424e1,7.94516e0,-1.36806e2}, + {-2.02617e0,-7.62769e-1,1.52768e1,-3.08728e1,-7.18354e1,7.05554e1}, + {3.5614e-1,-5.25061e-1,-6.86679e0,7.21247e0,4.25672e1,1.92463e0} + }}, + //Ho................................................... + CoeffMat{{ + {1.00004e0,-9.76747e-4,5.02266e-3,1.3853e-2,-4.89737e-2,-9.3934e-2}, + {-3.54744e-3,1.77664e0,1.04132e0,-1.22473e1,1.08741e1,4.94269e1}, + {2.21429e0,-3.09743e-1,-1.21254e1,3.59359e1,2.986e1,-1.11177e2}, + {-2.07893e0,-9.74603e-1,1.84979e1,-2.77795e1,-9.928e1,3.61173e1}, + {3.5809e-1,-5.08796e-1,-8.02233e0,5.91833e0,5.27374e1,1.53069e1} + }}, + //Er................................................... + CoeffMat{{ + {1.00007e0,-1.01117e-3,4.54349e-3,1.4172e-2,-4.56907e-2,-9.14834e-2}, + {-3.54805e-2,1.74848e0,1.74921e0,-1.20063e1,5.83106e0,4.40029e1}, + {2.3157e0,-3.56472e-2,-1.49988e1,3.37013e1,5.29727e1,-8.34879e1}, + {-2.12613e0,-1.19953e0,2.18208e1,-2.43067e1,-1.27958e2,-6.05797e-1}, + {3.57593e-1,-4.91966e-1,-9.2087e0,4.49274e0,6.33009e1,2.94631e1} + }}, + //Tm................................................... + CoeffMat{{ + {1.00009e0,-1.04219e-3,3.99979e-3,1.44301e-2,-4.18849e-2,-8.83308e-2}, + {-6.78832e-2,1.71323e0,2.49363e0,-1.16834e1,4.45978e-1,3.80313e1}, + {2.41433e0,2.59886e-1,-1.79765e1,3.11249e1,7.72617e1,-5.37152e1}, + {-2.16704e0,-1.4368e0,2.52369e1,-2.04429e1,-1.57828e2,-3.9614e1}, + {3.544e-1,-4.75067e-1,-1.04225e1,2.93269e0,7.42385e1,4.43866e1} + }}, + //Yb................................................... + CoeffMat{{ + {1.00012e0,-1.06927e-3,3.39091e-3,1.46208e-2,-3.75471e-2,-8.44506e-2}, + {-1.00611e-1,1.67068e0,3.27246e0,-1.12756e1,-5.27058e0,3.15152e1}, + {2.50948e0,5.76648e-1,-2.10477e1,2.81988e1,1.02665e2,-2.18989e1}, + {-2.20084e0,-1.68528e0,2.87325e1,-1.61828e1,-1.88802e2,-8.08394e1}, + {3.48226e-1,-4.58722e-1,-1.16589e1,1.23735e0,8.55166e1,6.00477e1} + }}, + //Lu................................................... + CoeffMat{{ + {1.00015e0,-1.09216e-3,2.71745e-3,1.47404e-2,-3.26779e-2,-7.98343e-2}, + {-1.33495e-1,1.62059e0,4.08307e0,-1.07793e1,-1.13042e1,2.44607e1}, + {2.60038e0,9.14519e-1,-2.42e1,2.49142e1,1.29108e2,1.19101e1}, + {-2.22661e0,-1.94391e0,3.22924e1,-1.15206e1,-2.20784e2,-1.24202e2}, + {3.38764e-1,-4.43546e-1,-1.29122e1,-5.9428e-1,9.70972e1,7.64128e1} + }}, + //Hf.................................................. + CoeffMat{{ + {1.00018e0,-1.11005e-3,1.98019e-3,1.47814e-2,-2.72814e-2,-7.44728e-2}, + {-1.66392e-1,1.56288e0,4.92359e0,-1.01915e1,-1.76472e1,1.68656e1}, + {2.68631e0,1.27279e0,-2.74231e1,2.12642e1,1.56537e2,4.76851e1}, + {-2.24352e0,-2.21096e0,3.59035e1,-6.45262e0,-2.53694e2,-1.69646e2}, + {3.25734e-1,-4.30391e-1,-1.41779e1,-2.56236e0,1.08949e2,9.34555e1} + }}, + //Ta.................................................. + CoeffMat{{ + {1.00021e0,-1.1227e-3,1.1875e-3,1.47467e-2,-2.14258e-2,-6.84575e-2}, + {-1.99145e-1,1.49739e0,5.79226e0,-9.50812e0,-2.42938e1,8.72312e0}, + {2.76656e0,1.65088e0,-3.07077e1,1.72391e1,1.84905e2,8.54178e1}, + {-2.25073e0,-2.48478e0,3.95541e1,-9.72452e-1,-2.87462e2,-2.17133e2}, + {3.0885e-1,-4.20088e-1,-1.54515e1,-4.66797e0,1.21043e2,1.11156e2} + }}, + //W................................................... + CoeffMat{{ + {1.00025e0,-1.12974e-3,3.25083e-4,1.46196e-2,-1.49818e-2,-6.15916e-2}, + {-2.31499e-1,1.42403e0,6.68301e0,-8.72894e0,-3.12005e1,8.13039e-2}, + {2.84004e0,2.04799e0,-3.40292e1,1.28426e1,2.14035e2,1.24907e2}, + {-2.24701e0,-2.76343e0,4.3216e1,4.91177e0,-3.2188e2,-2.66423e2}, + {2.87708e-1,-4.13579e-1,-1.67235e1,-6.90724e0,1.33305e2,1.29426e2} + }}, + //Re.................................................. + CoeffMat{{ + {1.00028e0,-1.13018e-3,-5.85208e-4,1.44008e-2,-8.12721e-3,-5.40793e-2}, + {-2.63328e-1,1.34285e0,7.59541e0,-7.85118e0,-3.83723e1,-9.07534e0}, + {2.90612e0,2.46267e0,-3.73829e1,8.06978e0,2.43918e2,1.66172e2}, + {-2.23164e0,-3.04426e0,4.68824e1,1.12007e1,-3.56918e2,-3.17506e2}, + {2.62068e-1,-4.12056e-1,-1.79911e1,-9.27922e0,1.45719e2,1.48255e2} + }}, + //Os.................................................. + CoeffMat{{ + {1.00032e0,-1.12442e-3,-1.55891e-3,1.40881e-2,-7.3226e-4,-4.57754e-2}, + {-2.94323e-1,1.25374e0,8.52167e0,-6.87354e0,-4.57537e1,-1.8691e1}, + {2.96357e0,2.89403e0,-4.07385e1,2.92128e0,2.74335e2,2.08986e2}, + {-2.20324e0,-3.32515e0,5.05197e1,1.78893e1,-3.92323e2,-3.70113e2}, + {2.31474e-1,-4.16536e-1,-1.92429e1,-1.1781e1,1.58197e2,1.67545e2} + }}, + //Ir................................................... + CoeffMat{{ + {1.00036e0,-1.11107e-3,-2.58076e-3,1.36692e-2,7.08911e-3,-3.6773e-2}, + {-3.24316e-1,1.15695e0,9.45947e0,-5.79675e0,-5.33328e1,-2.87528e1}, + {3.01159e0,3.33981e0,-4.4085e1,-2.59468e0,3.05217e2,2.53266e2}, + {-2.16093e0,-3.60244e0,5.41141e1,2.49638e1,-4.28002e2,-4.24126e2}, + {1.95645e-1,-4.28567e-1,-2.04742e1,-1.44066e1,1.70702e2,1.87249e2} + }}, + //Pt................................................... + CoeffMat{{ + {1.00039e0,-1.09011e-3,-3.65347e-3,1.31413e-2,1.53649e-2,-2.70335e-2}, + {-3.53028e-1,1.05257e0,1.04028e1,-4.62034e0,-6.10681e1,-3.92186e1}, + {3.04907e0,3.79831e0,-4.73979e1,-8.47391e0,3.36394e2,2.98829e2}, + {-2.10347e0,-3.87306e0,5.76385e1,3.24147e1,-4.63751e2,-4.79318e2}, + {1.54174e-1,-4.49506e-1,-2.16756e1,-1.71516e1,1.83161e2,2.07283e2} + }}, + //Au................................................... + CoeffMat{{ + {1.00043e0,-1.06114e-3,-4.7734e-3,1.25001e-2,2.40713e-2,-1.65728e-2}, + {-3.80207e-1,9.40778e-1,1.13464e1,-3.34423e0,-6.89259e1,-5.00549e1}, + {3.07493e0,4.26732e0,-5.06567e1,-1.47102e1,3.67721e2,3.45521e2}, + {-2.0297e0,-4.13334e0,6.10696e1,4.02298e1,-4.99399e2,-5.35495e2}, + {1.06694e-1,-4.80932e-1,-2.28393e1,-2.00106e1,1.95512e2,2.27574e2} + }}, + //Hg................................................. + CoeffMat{{ + {1.00047e0,-1.02377e-3,-5.93877e-3,1.17393e-2,3.32021e-2,-5.38025e-3}, + {-4.05577e-1,8.21834e-1,1.22846e1,-1.96909e0,-7.68675e1,-6.12207e1}, + {3.08808e0,4.74434e0,-5.38383e1,-2.12945e1,3.99035e2,3.93158e2}, + {-1.93842e0,-4.37923e0,6.43816e1,4.83938e1,-5.3475e2,-5.92428e2}, + {5.28167e-2,-5.24544e-1,-2.39565e1,-2.29768e1,2.07685e2,2.48037e2} + }}, + //Tl................................................. + CoeffMat{{ + {1.00051e0,-9.77691e-4,-7.14416e-3,1.08563e-2,4.27196e-2,6.50699e-3}, + {-4.28882e-1,6.96049e-1,1.32126e1,-4.94649e-1,-8.48638e1,-7.26907e1}, + {3.08745e0,5.22659e0,-5.69235e1,-2.82201e1,4.30208e2,4.41611e2}, + {-1.82849e0,-4.60634e0,6.7553e1,5.68935e1,-5.69648e2,-6.49946e2}, + {-7.81506e-3,-5.82177e-1,-2.50201e1,-2.60444e1,2.19621e2,2.68607e2} + }}, + //Pb................................................. + CoeffMat{{ + {1.00055e0,-9.22411e-4,-8.38967e-3,9.83942e-3,5.26394e-2,1.91419e-2}, + {-4.49833e-1,5.63865e-1,1.41246e1,1.0768e0,-9.28737e1,-8.44177e1}, + {3.07189e0,5.71077e0,-5.98888e1,-3.54715e1,4.61066e2,4.9067e2}, + {-1.69869e0,-4.80968e0,7.05579e1,6.57062e1,-6.03887e2,-7.07793e2}, + {-7.55851e-2,-6.55871e-1,-2.60213e1,-2.92042e1,2.31249e2,2.89188e2} + }}, + //Bi................................................. + CoeffMat{{ + {1.00059e0,-8.58152e-4,-9.6602e-3,8.70229e-3,6.28328e-2,3.23376e-2}, + {-4.68192e-1,4.25668e-1,1.50174e1,2.74763e0,-1.00883e2,-9.64012e1}, + {3.04043e0,6.19366e0,-6.27207e1,-4.30487e1,4.91532e2,5.40291e2}, + {-1.54799e0,-4.98431e0,7.3381e1,7.48252e1,-6.37365e2,-7.65885e2}, + {-1.5081e-1,-7.47659e-1,-2.69551e1,-3.24524e1,2.42528e2,3.09744e2} + }}, + //Po................................................. + CoeffMat{{ + {1.00062e0,-7.84582e-4,-1.09553e-2,7.43533e-3,7.33106e-2,4.61358e-2}, + {-4.83646e-1,2.8205e-1,1.58842e1,4.51493e0,-1.08846e2,-1.08589e2}, + {2.99185e0,6.67128e0,-6.53938e1,-5.09329e1,5.21419e2,5.90246e2}, + {-1.37512e0,-5.12439e0,7.59949e1,8.42238e1,-6.69862e2,-8.23945e2}, + {-2.33878e-1,-8.59879e-1,-2.78125e1,-3.57783e1,2.53381e2,3.30173e2} + }}, + //At................................................. + CoeffMat{{ + {1.00066e0,-7.01484e-4,-1.22684e-2,6.03454e-3,8.40298e-2,6.05018e-2}, + {-4.95899e-1,1.3369e-1,1.67195e1,6.37574e0,-1.16723e2,-1.20936e2}, + {2.92499e0,7.13921e0,-6.78861e1,-5.91051e1,5.50562e2,6.40337e2}, + {-1.17889e0,-5.22361e0,7.8376e1,9.38748e1,-7.01185e2,-8.81728e2}, + {-3.25146e-1,-9.95032e-1,-2.85858e1,-3.91713e1,2.63738e2,3.50385e2} + }}, + //Rn................................................. + CoeffMat{{ + {1.0007e0,-6.08729e-4,-1.35873e-2,4.50061e-3,9.48991e-2,7.53277e-2}, + {-5.04647e-1,-1.86445e-2,1.75175e1,8.32729e0,-1.24477e2,-1.33405e2}, + {2.83871e0,7.59269e0,-7.01766e1,-6.75463e1,5.78809e2,6.90388e2}, + {-9.58143e-1,-5.27523e0,8.0502e1,1.03751e2,-7.31155e2,-9.39013e2}, + {-4.24964e-1,-1.15578e0,-2.92679e1,-4.26205e1,2.73537e2,3.70298e2} + }}, + //Fr................................................. + CoeffMat{{ + {1.00073e0,-5.06594e-4,-1.49195e-2,2.82752e-3,1.05984e-1,9.07179e-2}, + {-5.09498e-1,-1.74114e-1,1.82691e1,1.0364e1,-1.32042e2,-1.45912e2}, + {2.73154e0,8.02658e0,-7.22315e1,-7.62278e1,6.05902e2,7.40063e2}, + {-7.11399e-1,-5.27209e0,8.23375e1,1.13815e2,-7.59479e2,-9.95409e2}, + {-5.33773e-1,-1.34491e0,-2.98474e1,-4.61115e1,2.82675e2,3.89772e2} + }}, + //Ra................................................. + CoeffMat{{ + {1.00076e0,-3.94938e-4,-1.62461e-2,1.01616e-3,1.17156e-1,1.06529e-1}, + {-5.10179e-1,-3.3169e-1,1.89693e1,1.24811e1,-1.39382e2,-1.58412e2}, + {2.60247e0,8.43503e0,-7.4033e1,-8.51228e1,6.31692e2,7.89166e2}, + {-4.37673e-1,-5.2062e0,8.38642e1,1.2403e2,-7.85984e2,-1.05067e3}, + {-6.5185e-1,-1.56551e0,-3.03186e1,-4.96307e1,2.91091e2,4.08719e2} + }}, + //Ac................................................. + CoeffMat{{ + {1.00079e0,-2.74014e-4,-1.75614e-2,-9.33046e-4,1.28369e-1,1.22714e-1}, + {-5.06317e-1,-4.90282e-1,1.96102e1,1.46725e1,-1.4644e2,-1.70838e2}, + {2.45014e0,8.81187e0,-7.5553e1,-9.42005e1,6.55963e2,8.37419e2}, + {-1.35635e-1,-5.06921e0,8.50532e1,1.34356e2,-8.10425e2,-1.10448e3}, + {-7.79579e-1,-1.82077e0,-3.06725e1,-5.31628e1,2.987e2,4.2702e2} + }}, + //Th................................................ + CoeffMat{{ + {1.00082e0,-1.43895e-4,-1.8859e-2,-3.02404e-3,1.39581e-1,1.39236e-1}, + {-4.97537e-1,-6.48585e-1,2.01838e1,1.69296e1,-1.53157e2,-1.83108e2}, + {2.27321e0,9.15011e0,-7.6763e1,-1.0342e2,6.78484e2,8.84496e2}, + {1.9599e-1,-4.85189e0,8.5876e1,1.44742e2,-8.32543e2,-1.15645e3}, + {-9.17319e-1,-2.1142e0,-3.09003e1,-5.66889e1,3.05413e2,4.44542e2} + }}, + //Pa................................................. + CoeffMat{{ + {1.00085e0,-5.30411e-6,-2.01239e-2,-5.2478e-3,1.50683e-1,1.55963e-1}, + {-4.83526e-1,-8.05262e-1,2.06846e1,1.92461e1,-1.59496e2,-1.95177e2}, + {2.07059e0,9.44256e0,-7.76442e1,-1.12749e2,6.99103e2,9.30193e2}, + {5.58224e-1,-4.54471e0,8.63141e1,1.55146e2,-8.52163e2,-1.20633e3}, + {-1.06534e0,-2.44938e0,-3.09967e1,-6.01932e1,3.11168e2,4.61192e2} + }}, + //U.................................................... + CoeffMat{{ + {1.00087e0,1.42049e-4,-2.13455e-2,-7.61031e-3,1.61591e-1,1.72814e-1}, + {-4.6392e-1,-9.58725e-1,2.11056e1,2.16134e1,-1.65402e2,-2.06974e2}, + {1.841e0,9.68106e0,-7.81727e1,-1.22145e2,7.17614e2,9.74218e2}, + {9.52222e-1,-4.13715e0,8.63446e1,1.65515e2,-8.69054e2,-1.25378e3}, + {-1.22395e0,-2.83024e0,-3.09548e1,-6.36564e1,3.15887e2,4.76849e2} + }}}; + // clang-format on + static_assert(std::size(mott_coeffs) == MottElementData::num_mott_elements, + "wrong number of Mott coefficient elements"); + + int index = z.unchecked_get() - 1; + CELER_VALIDATE( + index >= 0 && index < int{MottElementData::num_mott_elements}, + << "atomic number " << z.get() + << " is out of range for Coulomb scattering model Mott coefficients " + "(must be less than " + << MottElementData::num_mott_elements << ")"); + + return mott_coeffs[index]; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/params/WentzelOKVIParams.hh b/src/celeritas/em/params/WentzelOKVIParams.hh new file mode 100644 index 0000000000..14b93495de --- /dev/null +++ b/src/celeritas/em/params/WentzelOKVIParams.hh @@ -0,0 +1,85 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 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/params/WentzelOKVIParams.hh +//---------------------------------------------------------------------------// +#pragma once + +#include + +#include "corecel/data/CollectionMirror.hh" +#include "corecel/data/ParamsDataInterface.hh" +#include "celeritas/em/data/WentzelOKVIData.hh" +#include "celeritas/mat/IsotopeView.hh" +#include "celeritas/phys/AtomicNumber.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +class MaterialParams; +struct ImportData; + +//---------------------------------------------------------------------------// +/*! + * Construct and store shared Coulomb and multiple scattering data. + * + * This data is used by both the single Coulomb scattering and Wentzel VI + * multiple scattering models. + */ +class WentzelOKVIParams final : public ParamsDataInterface +{ + public: + //!@{ + //! \name Type aliases + using SPConstMaterials = std::shared_ptr; + //!@} + + public: + struct Options + { + //! Use combined single and multiple scattering + bool is_combined{true}; + //! Polar angle limit between single and multiple scattering + real_type polar_angle_limit{constants::pi}; + //! Factor for dynamic computation of angular limit between SS and MSC + real_type angle_limit_factor{1}; + //! User defined screening factor + real_type screening_factor{1}; + //! Nuclear form factor model + NuclearFormFactorType form_factor{NuclearFormFactorType::exponential}; + }; + + public: + // Construct if Wentzel VI or Coulomb is present, else return nullptr + static std::shared_ptr + from_import(ImportData const& data, SPConstMaterials materials); + + // Construct from material data and options + WentzelOKVIParams(SPConstMaterials materials, Options options); + + //! Access Wentzel OK&VI data on the host + HostRef const& host_ref() const final { return data_.host_ref(); } + + //! Access Wentzel OK&VI data on the device + DeviceRef const& device_ref() const final { return data_.device_ref(); } + + private: + // Host/device storage and reference + CollectionMirror data_; + + // Construct per-element data (loads Mott coefficients) + void build_data(HostVal& host_data, + MaterialParams const& materials); + + // Retrieve matrix of interpolated Mott coefficients + static MottElementData::MottCoeffMatrix + get_mott_coeff_matrix(AtomicNumber z); + + // Calculate the nuclear form prefactor + static real_type calc_nuclear_form_prefactor(IsotopeView const& iso); +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/params/WentzelVIMscParams.cc b/src/celeritas/em/params/WentzelVIMscParams.cc index 6f7048a994..090aae4189 100644 --- a/src/celeritas/em/params/WentzelVIMscParams.cc +++ b/src/celeritas/em/params/WentzelVIMscParams.cc @@ -29,10 +29,7 @@ WentzelVIMscParams::from_import(ParticleParams const& particles, MaterialParams const& materials, ImportData const& data) { - auto is_wentzel = [](ImportMscModel const& imm) { - return imm.model_class == ImportModelClass::wentzel_vi_uni; - }; - if (!std::any_of(data.msc_models.begin(), data.msc_models.end(), is_wentzel)) + if (!has_msc_model(data, ImportModelClass::wentzel_vi_uni)) { // No WentzelVI MSC present return nullptr; diff --git a/src/celeritas/em/params/detail/MscParamsHelper.cc b/src/celeritas/em/params/detail/MscParamsHelper.cc index 5b80c43988..52ea078907 100644 --- a/src/celeritas/em/params/detail/MscParamsHelper.cc +++ b/src/celeritas/em/params/detail/MscParamsHelper.cc @@ -43,7 +43,7 @@ MscParamsHelper::MscParamsHelper(ParticleParams const& particles, /*! * Validate and save MSC IDs. */ -void MscParamsHelper::build_ids(MscIds* ids) const +void MscParamsHelper::build_ids(CoulombIds* ids) const { ids->electron = particles_.find(pdg::electron()); ids->positron = particles_.find(pdg::positron()); diff --git a/src/celeritas/em/params/detail/MscParamsHelper.hh b/src/celeritas/em/params/detail/MscParamsHelper.hh index 978c1fe8bd..0059894f33 100644 --- a/src/celeritas/em/params/detail/MscParamsHelper.hh +++ b/src/celeritas/em/params/detail/MscParamsHelper.hh @@ -11,7 +11,7 @@ #include "corecel/Types.hh" #include "corecel/data/Collection.hh" -#include "celeritas/em/data/MscData.hh" +#include "celeritas/em/data/CommonCoulombData.hh" #include "celeritas/io/ImportModel.hh" namespace celeritas @@ -42,7 +42,7 @@ class MscParamsHelper VecImportMscModel const&, ImportModelClass); - void build_ids(MscIds* ids) const; + void build_ids(CoulombIds* ids) const; void build_xs(XsValues*, Values*) const; private: diff --git a/src/celeritas/em/process/CoulombScatteringProcess.cc b/src/celeritas/em/process/CoulombScatteringProcess.cc index 5eac088a7e..aafc0f2ee3 100644 --- a/src/celeritas/em/process/CoulombScatteringProcess.cc +++ b/src/celeritas/em/process/CoulombScatteringProcess.cc @@ -19,13 +19,10 @@ namespace celeritas /*! * Construct from host data. */ -CoulombScatteringProcess::CoulombScatteringProcess( - SPConstParticles particles, - SPConstMaterials materials, - SPConstImported process_data, - CoulombScatteringModel::Options const& options) +CoulombScatteringProcess::CoulombScatteringProcess(SPConstParticles particles, + SPConstImported process_data, + Options const& options) : particles_(std::move(particles)) - , materials_(std::move(materials)) , imported_(process_data, particles_, ImportProcessClass::coulomb_scat, @@ -33,8 +30,6 @@ CoulombScatteringProcess::CoulombScatteringProcess( , options_(options) { CELER_EXPECT(particles_); - CELER_EXPECT(materials_); - CELER_EXPECT(options_); } //---------------------------------------------------------------------------// @@ -45,7 +40,7 @@ auto CoulombScatteringProcess::build_models(ActionIdIter start_id) const -> VecModel { return {std::make_shared( - *start_id++, *particles_, *materials_, options_, imported_.processes())}; + *start_id++, *particles_, imported_.processes())}; } //---------------------------------------------------------------------------// diff --git a/src/celeritas/em/process/CoulombScatteringProcess.hh b/src/celeritas/em/process/CoulombScatteringProcess.hh index a2f8db943d..92038b2ab4 100644 --- a/src/celeritas/em/process/CoulombScatteringProcess.hh +++ b/src/celeritas/em/process/CoulombScatteringProcess.hh @@ -12,7 +12,6 @@ #include "celeritas/em/data/CoulombScatteringData.hh" #include "celeritas/em/model/CoulombScatteringModel.hh" #include "celeritas/io/ImportParameters.hh" -#include "celeritas/mat/MaterialParams.hh" #include "celeritas/phys/Applicability.hh" #include "celeritas/phys/ImportedProcessAdapter.hh" #include "celeritas/phys/Process.hh" @@ -30,16 +29,20 @@ class CoulombScatteringProcess : public Process //!@{ //! \name Type aliases using SPConstParticles = std::shared_ptr; - using SPConstMaterials = std::shared_ptr; using SPConstImported = std::shared_ptr; //!@} + struct Options + { + //! Whether to use integral method to sample interaction length + bool use_integral_xs{true}; + }; + public: //! Construct from Coulomb scattering data CoulombScatteringProcess(SPConstParticles particles, - SPConstMaterials materials, SPConstImported process_data, - CoulombScatteringModel::Options const& options); + Options const& options); //! Construct the models associated with this process VecModel build_models(ActionIdIter start_id) const final; @@ -55,9 +58,8 @@ class CoulombScatteringProcess : public Process private: SPConstParticles particles_; - SPConstMaterials materials_; ImportedProcessAdapter imported_; - CoulombScatteringModel::Options options_; + Options options_; }; //---------------------------------------------------------------------------// diff --git a/src/celeritas/em/xs/MottRatioCalculator.hh b/src/celeritas/em/xs/MottRatioCalculator.hh index 28d626a639..68f9575aae 100644 --- a/src/celeritas/em/xs/MottRatioCalculator.hh +++ b/src/celeritas/em/xs/MottRatioCalculator.hh @@ -10,7 +10,7 @@ #include "corecel/Macros.hh" #include "corecel/Types.hh" #include "corecel/math/ArrayUtils.hh" -#include "celeritas/em/data/CoulombScatteringData.hh" +#include "celeritas/em/data/WentzelOKVIData.hh" #include "celeritas/grid/PolyEvaluator.hh" namespace celeritas @@ -31,14 +31,13 @@ class MottRatioCalculator public: //! Construct with state data inline CELER_FUNCTION - MottRatioCalculator(CoulombScatteringElementData const& element_data, - real_type beta); + MottRatioCalculator(MottElementData const& element_data, real_type beta); //! Ratio of Mott and Rutherford cross sections inline CELER_FUNCTION real_type operator()(real_type cos_t) const; private: - CoulombScatteringElementData const& element_data_; + MottElementData const& element_data_; real_type beta_; }; @@ -49,8 +48,8 @@ class MottRatioCalculator * Construct with state data. */ CELER_FUNCTION -MottRatioCalculator::MottRatioCalculator( - CoulombScatteringElementData const& element_data, real_type beta) +MottRatioCalculator::MottRatioCalculator(MottElementData const& element_data, + real_type beta) : element_data_(element_data), beta_(beta) { CELER_EXPECT(0 <= beta_ && beta_ < 1); @@ -80,12 +79,14 @@ real_type MottRatioCalculator::operator()(real_type cos_theta) const real_type beta0 = beta_ - beta_shift; // Evaluate polynomial of powers of beta0 and fcos_t - CoulombScatteringElementData::ThetaArray theta_coeffs; + MottElementData::ThetaArray theta_coeffs; for (auto i : range(theta_coeffs.size())) { theta_coeffs[i] = PolyEvaluator(element_data_.mott_coeff[i])(beta0); } - return PolyEvaluator(theta_coeffs)(fcos_t); + real_type result = PolyEvaluator(theta_coeffs)(fcos_t); + CELER_ENSURE(result >= 0); + return result; } //---------------------------------------------------------------------------// diff --git a/src/celeritas/em/xs/WentzelHelper.hh b/src/celeritas/em/xs/WentzelHelper.hh index be9d92155c..e9a6679f0a 100644 --- a/src/celeritas/em/xs/WentzelHelper.hh +++ b/src/celeritas/em/xs/WentzelHelper.hh @@ -10,6 +10,9 @@ #include "corecel/Macros.hh" #include "corecel/Types.hh" #include "corecel/math/Algorithms.hh" +#include "celeritas/em/data/CommonCoulombData.hh" +#include "celeritas/em/data/WentzelOKVIData.hh" +#include "celeritas/mat/MaterialView.hh" #include "celeritas/phys/AtomicNumber.hh" #include "celeritas/phys/ParticleTrackView.hh" @@ -31,17 +34,21 @@ class WentzelHelper public: //!@{ //! \name Type aliases - using MomentumSq = units::MevMomentumSq; + using Charge = units::ElementaryCharge; using Energy = units::MevEnergy; using Mass = units::MevMass; + using MomentumSq = units::MevMomentumSq; //!@} public: // Construct from particle and material properties - inline CELER_FUNCTION WentzelHelper(ParticleTrackView const& particle, - AtomicNumber target_z, - CoulombScatteringRef const& data, - Energy cutoff); + inline CELER_FUNCTION + WentzelHelper(ParticleTrackView const& particle, + MaterialView const& material, + AtomicNumber target_z, + NativeCRef const& wentzel, + CoulombIds const& ids, + Energy cutoff); //! Get the target atomic number CELER_FUNCTION AtomicNumber atomic_number() const { return target_z_; } @@ -52,26 +59,45 @@ class WentzelHelper return screening_coefficient_; } + //! Get the Mott factor + CELER_FUNCTION real_type mott_factor() const { return mott_factor_; } + + //! Get the multiplicative factor for the cross section + CELER_FUNCTION real_type kin_factor() const { return kin_factor_; } + //! Get the maximum scattering angle off of electrons - CELER_FUNCTION real_type costheta_max_electron() const + CELER_FUNCTION real_type cos_thetamax_electron() const + { + return cos_thetamax_elec_; + } + + //! Get the maximum scattering angle off of a nucleus + CELER_FUNCTION real_type cos_thetamax_nuclear() const { - return cos_t_max_elec_; + return cos_thetamax_nuc_; } // The ratio of electron to total cross section for Coulomb scattering - inline CELER_FUNCTION real_type calc_xs_ratio() const; + inline CELER_FUNCTION real_type calc_xs_ratio(real_type cos_thetamin, + real_type cos_thetamax) const; + + // Calculate the electron cross section for Coulomb scattering + inline CELER_FUNCTION real_type + calc_xs_electron(real_type cos_thetamin, real_type cos_thetamax) const; + + // Calculate the nuclear cross section for Coulomb scattering + inline CELER_FUNCTION real_type + calc_xs_nuclear(real_type cos_thetamin, real_type cos_thetamax) const; private: //// DATA //// - // Target atomic number AtomicNumber const target_z_; - - // Moliere screening coefficient real_type screening_coefficient_; - - // Cosine of the maximum scattering angle off of electrons - real_type cos_t_max_elec_; + real_type kin_factor_; + real_type mott_factor_; + real_type cos_thetamax_elec_; + real_type cos_thetamax_nuc_; //// HELPER FUNCTIONS //// @@ -82,9 +108,23 @@ class WentzelHelper // Calculate the screening coefficient R^2 for electrons CELER_CONSTEXPR_FUNCTION real_type screen_r_sq_elec() const; + // Calculate the multiplicative factor for the cross section + inline CELER_FUNCTION real_type + calc_kin_factor(ParticleTrackView const&) const; + // Calculate the (cosine of) the maximum scattering angle off of electrons - inline CELER_FUNCTION real_type calc_costheta_max_electron( - ParticleTrackView const&, CoulombScatteringRef const&, Energy) const; + inline CELER_FUNCTION real_type calc_cos_thetamax_electron( + ParticleTrackView const&, CoulombIds const&, Energy) const; + + // Calculate the (cosine of) the maximum scattering angle off of a nucleus + inline CELER_FUNCTION real_type calc_cos_thetamax_nuclear( + ParticleTrackView const&, + MaterialView const& material, + NativeCRef const& wentzel) const; + + // Calculate the common factor in the electron and nuclear cross section + inline CELER_FUNCTION real_type calc_xs_factor(real_type cos_thetamin, + real_type cos_thetamax) const; }; //---------------------------------------------------------------------------// @@ -95,31 +135,76 @@ class WentzelHelper */ CELER_FUNCTION WentzelHelper::WentzelHelper(ParticleTrackView const& particle, + MaterialView const& material, AtomicNumber target_z, - CoulombScatteringRef const& data, + NativeCRef const& wentzel, + CoulombIds const& ids, Energy cutoff) : target_z_(target_z) , screening_coefficient_(this->calc_screening_coefficient(particle) - * data.screening_factor) - , cos_t_max_elec_(this->calc_costheta_max_electron(particle, data, cutoff)) + * wentzel.params.screening_factor) + , kin_factor_(this->calc_kin_factor(particle)) + , mott_factor_(particle.particle_id() == ids.electron + ? 1 + real_type(2e-4) * ipow<2>(target_z_.get()) + : 1) + , cos_thetamax_elec_( + this->calc_cos_thetamax_electron(particle, ids, cutoff)) + , cos_thetamax_nuc_( + this->calc_cos_thetamax_nuclear(particle, material, wentzel)) { CELER_EXPECT(screening_coefficient_ > 0); - CELER_EXPECT(cos_t_max_elec_ >= -1 && cos_t_max_elec_ <= 1); + CELER_EXPECT(cos_thetamax_elec_ >= -1 && cos_thetamax_elec_ <= 1); + CELER_EXPECT(cos_thetamax_nuc_ >= -1 && cos_thetamax_nuc_ <= 1); } //---------------------------------------------------------------------------// /*! * Ratio of electron cross section to total (nuclear + electron) cross section. */ -CELER_FUNCTION real_type WentzelHelper::calc_xs_ratio() const +CELER_FUNCTION real_type WentzelHelper::calc_xs_ratio( + real_type cos_thetamin, real_type cos_thetamax) const +{ + real_type xs_elec = this->calc_xs_electron(cos_thetamin, cos_thetamax); + return xs_elec + / (xs_elec + this->calc_xs_nuclear(cos_thetamin, cos_thetamax)); +} + +//---------------------------------------------------------------------------// +/*! + * Calculate the electron cross section for Coulomb scattering. + */ +CELER_FUNCTION real_type WentzelHelper::calc_xs_electron( + real_type cos_thetamin, real_type cos_thetamax) const +{ + cos_thetamin = max(cos_thetamin, cos_thetamax_elec_); + cos_thetamax = max(cos_thetamax, cos_thetamax_elec_); + if (cos_thetamin <= cos_thetamax) + { + return 0; + } + return this->calc_xs_factor(cos_thetamin, cos_thetamax); +} + +//---------------------------------------------------------------------------// +/*! + * Calculate the nuclear cross section for Coulomb scattering. + */ +CELER_FUNCTION real_type WentzelHelper::calc_xs_nuclear( + real_type cos_thetamin, real_type cos_thetamax) const { - // Calculating only reduced cross sections by elimination mutual factors - // in the ratio. - real_type nuc_xsec = target_z_.get() / (1 + screening_coefficient_); - real_type elec_xsec = (1 - cos_t_max_elec_) - / (1 - cos_t_max_elec_ + 2 * screening_coefficient_); + return target_z_.get() * this->calc_xs_factor(cos_thetamin, cos_thetamax); +} - return elec_xsec / (nuc_xsec + elec_xsec); +//---------------------------------------------------------------------------// +/*! + * Calculate the common factor in the electron and nuclear cross section. + */ +CELER_FUNCTION real_type WentzelHelper::calc_xs_factor( + real_type cos_thetamin, real_type cos_thetamax) const +{ + return kin_factor_ * mott_factor_ * (cos_thetamin - cos_thetamax) + / ((1 - cos_thetamin + 2 * screening_coefficient_) + * (1 - cos_thetamax + 2 * screening_coefficient_)); } //---------------------------------------------------------------------------// @@ -178,6 +263,32 @@ CELER_CONSTEXPR_FUNCTION real_type WentzelHelper::screen_r_sq_elec() const .value(); } +//---------------------------------------------------------------------------// +/*! + * Calculate the multiplicative factor for the cross section. + * + * This calculates the factor + * \f[ + f = \frac{2 \pi m_e^2 r_e^2 Z q^2}{\beta^2 p^2}, + * \f] + * where \f$ m_e, r_e, Z, q, \beta \f$, and \f$ p \f$ are the electron mass, + * classical electron radius, atomic number of the target atom, charge, + * relativistic speed, and momentum of the incident particle, respectively. + */ +CELER_FUNCTION real_type +WentzelHelper::calc_kin_factor(ParticleTrackView const& particle) const +{ + real_type constexpr twopi_mrsq + = 2 * constants::pi + * ipow<2>(native_value_to(constants::electron_mass).value() + * constants::r_electron); + + return twopi_mrsq * target_z_.get() + * ipow<2>(value_as(particle.charge())) + / (particle.beta_sq() + * value_as(particle.momentum_sq())); +} + //---------------------------------------------------------------------------// /*! * Calculate the maximum scattering angle off the target's electrons. @@ -185,15 +296,13 @@ CELER_CONSTEXPR_FUNCTION real_type WentzelHelper::screen_r_sq_elec() const * This calculates the cosine of the maximum polar angle that the incident * particle can scatter off of the target's electrons. */ -CELER_FUNCTION real_type -WentzelHelper::calc_costheta_max_electron(ParticleTrackView const& particle, - CoulombScatteringRef const& data, - Energy cutoff) const +CELER_FUNCTION real_type WentzelHelper::calc_cos_thetamax_electron( + ParticleTrackView const& particle, CoulombIds const& ids, Energy cutoff) const { real_type inc_energy = value_as(particle.energy()); real_type mass = value_as(particle.mass()); - real_type max_energy = particle.particle_id() == data.ids.electron + real_type max_energy = particle.particle_id() == ids.electron ? real_type{0.5} * inc_energy : inc_energy; real_type final_energy = inc_energy @@ -210,5 +319,26 @@ WentzelHelper::calc_costheta_max_electron(ParticleTrackView const& particle, return 0; } +//---------------------------------------------------------------------------// +/*! + * Calculate the maximum scattering angle off the target nucleus. + */ +CELER_FUNCTION real_type WentzelHelper::calc_cos_thetamax_nuclear( + ParticleTrackView const& particle, + MaterialView const& material, + NativeCRef const& wentzel) const +{ + if (wentzel.params.is_combined) + { + CELER_ASSERT(material.material_id() < wentzel.inv_mass_cbrt_sq.size()); + return max(wentzel.params.costheta_limit, + 1 + - wentzel.params.a_sq_factor + * wentzel.inv_mass_cbrt_sq[material.material_id()] + / value_as(particle.momentum_sq())); + } + return wentzel.params.costheta_limit; +} + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/em/xs/WentzelTransportXsCalculator.hh b/src/celeritas/em/xs/WentzelTransportXsCalculator.hh index fa8e92b0b8..5bef790a91 100644 --- a/src/celeritas/em/xs/WentzelTransportXsCalculator.hh +++ b/src/celeritas/em/xs/WentzelTransportXsCalculator.hh @@ -30,7 +30,6 @@ class WentzelTransportXsCalculator //!@{ //! \name Type aliases using XsUnits = units::Native; // [len^2] - using Charge = units::ElementaryCharge; using Mass = units::MevMass; using MomentumSq = units::MevMomentumSq; //!@} @@ -42,26 +41,23 @@ class WentzelTransportXsCalculator WentzelHelper const& helper); // Calculate the transport cross section for the given angle [len^2] - inline CELER_FUNCTION real_type operator()(real_type costheta_max) const; + inline CELER_FUNCTION real_type operator()(real_type cos_thetamax) const; private: //// DATA //// AtomicNumber z_; real_type screening_coeff_; - real_type costheta_max_elec_; + real_type cos_thetamax_elec_; real_type beta_sq_; - real_type xs_factor_; + real_type kin_factor_; //// HELPER FUNCTIONS //// - // Calculate the multiplicative factor for the transport cross section - real_type calc_xs_factor(ParticleTrackView const&) const; - // Calculate xs contribution from scattering off electrons or nucleus - real_type calc_xs_contribution(real_type costheta_max) const; + real_type calc_xs_contribution(real_type cos_thetamax) const; - //! Limit on (1 - \c costheta_max) / \c screening_coeff + //! Limit on (1 - \c cos_thetamax) / \c screening_coeff static CELER_CONSTEXPR_FUNCTION real_type limit() { return 0.1; } }; @@ -72,7 +68,7 @@ class WentzelTransportXsCalculator * Construct with particle and precalculatad Wentzel data. * * \c beta_sq should be calculated from the incident particle energy and mass. - * \c screening_coeff and \c costheta_max_elec are calculated using the Wentzel + * \c screening_coeff and \c cos_thetamax_elec are calculated using the Wentzel * OK and VI model in \c WentzelHelper and depend on properties of the incident * particle, the energy cutoff in the current material, and the target element. */ @@ -81,9 +77,9 @@ WentzelTransportXsCalculator::WentzelTransportXsCalculator( ParticleTrackView const& particle, WentzelHelper const& helper) : z_(helper.atomic_number()) , screening_coeff_(2 * helper.screening_coefficient()) - , costheta_max_elec_(helper.costheta_max_electron()) + , cos_thetamax_elec_(helper.cos_thetamax_electron()) , beta_sq_(particle.beta_sq()) - , xs_factor_(this->calc_xs_factor(particle)) + , kin_factor_(helper.kin_factor()) { } @@ -92,56 +88,31 @@ WentzelTransportXsCalculator::WentzelTransportXsCalculator( * Calculate the transport cross section for the given angle [len^2]. */ CELER_FUNCTION real_type -WentzelTransportXsCalculator::operator()(real_type costheta_max) const +WentzelTransportXsCalculator::operator()(real_type cos_thetamax) const { - CELER_EXPECT(costheta_max <= 1); + CELER_EXPECT(cos_thetamax <= 1); // Sum xs contributions from scattering off electrons and nucleus - real_type xs_nuc = this->calc_xs_contribution(costheta_max); - real_type xs_elec = costheta_max_elec_ > costheta_max - ? this->calc_xs_contribution(costheta_max_elec_) + real_type xs_nuc = this->calc_xs_contribution(cos_thetamax); + real_type xs_elec = cos_thetamax_elec_ > cos_thetamax + ? this->calc_xs_contribution(cos_thetamax_elec_) : xs_nuc; - real_type result = xs_factor_ * (xs_elec + z_.get() * xs_nuc); + real_type result = kin_factor_ * (xs_elec + z_.get() * xs_nuc); CELER_ENSURE(result >= 0); return result; } -//---------------------------------------------------------------------------// -/*! - * Calculate the multiplicative factor for the transport cross section. - * - * This calculates the factor - * \f[ - f = \frac{2 \pi m_e^2 r_e^2 Z q^2}{\beta^2 p^2}, - * \f] - * where \f$ m_e, r_e, Z, q, \beta \f$, and \f$ p \f$ are the electron mass, - * classical electron radius, atomic number of the target atom, charge, - * relativistic speed, and momentum of the incident particle, respectively. - */ -CELER_FUNCTION real_type WentzelTransportXsCalculator::calc_xs_factor( - ParticleTrackView const& particle) const -{ - real_type constexpr twopi_mrsq - = 2 * constants::pi - * ipow<2>(native_value_to(constants::electron_mass).value() - * constants::r_electron); - - return twopi_mrsq * z_.get() * ipow<2>(value_as(particle.charge())) - / (particle.beta_sq() - * value_as(particle.momentum_sq())); -} - //---------------------------------------------------------------------------// /*! * Calculate contribution to xs from scattering off electrons or nucleus. */ CELER_FUNCTION real_type WentzelTransportXsCalculator::calc_xs_contribution( - real_type costheta_max) const + real_type cos_thetamax) const { real_type result; real_type const spin = real_type(0.5); - real_type x = (1 - costheta_max) / screening_coeff_; + real_type x = (1 - cos_thetamax) / screening_coeff_; if (x < WentzelTransportXsCalculator::limit()) { real_type x_sq = ipow<2>(x); diff --git a/src/celeritas/ext/GeantImporter.cc b/src/celeritas/ext/GeantImporter.cc index d28bf1a4a1..a7b9f6e445 100644 --- a/src/celeritas/ext/GeantImporter.cc +++ b/src/celeritas/ext/GeantImporter.cc @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include @@ -327,6 +328,28 @@ to_msc_step_algorithm(G4MscStepLimitType const& msc_step_algorithm) CELER_ASSERT_UNREACHABLE(); } +//---------------------------------------------------------------------------// +/*! + * Safely switch from G4NuclearFormfactorType [G4NuclearFormfactorType.hh] to + * NuclearFormFactorType. + */ +NuclearFormFactorType +to_form_factor_type(G4NuclearFormfactorType const& form_factor_type) +{ + switch (form_factor_type) + { + case G4NuclearFormfactorType::fNoneNF: + return NuclearFormFactorType::none; + case G4NuclearFormfactorType::fExponentialNF: + return NuclearFormFactorType::exponential; + case G4NuclearFormfactorType::fGaussianNF: + return NuclearFormFactorType::gaussian; + case G4NuclearFormfactorType::fFlatNF: + return NuclearFormFactorType::flat; + } + CELER_ASSERT_UNREACHABLE(); +} + //---------------------------------------------------------------------------// /*! * Return a populated \c ImportParticle vector. @@ -919,8 +942,11 @@ ImportEmParameters import_em_parameters() #else CELER_DISCARD(len_scale); #endif + import.msc_theta_limit = g4.MscThetaLimit(); + import.angle_limit_factor = g4.FactorForAngleLimit(); import.apply_cuts = g4.ApplyCuts(); import.screening_factor = g4.ScreeningFactor(); + import.form_factor = to_form_factor_type(g4.NuclearFormfactorType()); CELER_ENSURE(import); return import; diff --git a/src/celeritas/ext/GeantPhysicsOptions.hh b/src/celeritas/ext/GeantPhysicsOptions.hh index 30a1fb0f87..1160cb02b7 100644 --- a/src/celeritas/ext/GeantPhysicsOptions.hh +++ b/src/celeritas/ext/GeantPhysicsOptions.hh @@ -8,6 +8,7 @@ #pragma once #include "corecel/Types.hh" +#include "celeritas/Constants.hh" #include "celeritas/Quantities.hh" namespace celeritas @@ -121,8 +122,14 @@ struct GeantPhysicsOptions double msc_safety_factor{0.6}; //! Lambda limit for MSC models [len] double msc_lambda_limit{0.1 * units::centimeter}; + //! Polar angle limii between single and multiple Coulomb scattering + double msc_theta_limit{constants::pi}; + //! Factor for dynamic computation of angular limit between SS and MSC + double angle_limit_factor{1}; //! Step limit algorithm for MSC models MscStepLimitAlgorithm msc_step_algorithm{MscStepLimitAlgorithm::safety}; + //! Nuclear form factor model for Coulomm scattering + NuclearFormFactorType form_factor{NuclearFormFactorType::exponential}; //!@} //! Print detailed Geant4 output @@ -158,7 +165,10 @@ operator==(GeantPhysicsOptions const& a, GeantPhysicsOptions const& b) && 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.msc_theta_limit == b.msc_theta_limit + && a.angle_limit_factor == b.angle_limit_factor && a.msc_step_algorithm == b.msc_step_algorithm + && a.form_factor == b.form_factor && a.verbose == b.verbose; // clang-format on } diff --git a/src/celeritas/ext/GeantPhysicsOptionsIO.json.cc b/src/celeritas/ext/GeantPhysicsOptionsIO.json.cc index 01002d2040..bc7cd6acbe 100644 --- a/src/celeritas/ext/GeantPhysicsOptionsIO.json.cc +++ b/src/celeritas/ext/GeantPhysicsOptionsIO.json.cc @@ -70,6 +70,19 @@ void to_json(nlohmann::json& j, MscStepLimitAlgorithm const& value) j = std::string{to_cstring(value)}; } +void from_json(nlohmann::json const& j, NuclearFormFactorType& value) +{ + static auto const from_string + = StringEnumMapper::from_cstring_func( + to_cstring, "form factor"); + value = from_string(j.get()); +} + +void to_json(nlohmann::json& j, NuclearFormFactorType const& value) +{ + j = std::string{to_cstring(value)}; +} + //---------------------------------------------------------------------------// /*! * Read options from JSON. @@ -108,7 +121,10 @@ void from_json(nlohmann::json const& j, GeantPhysicsOptions& options) GPO_LOAD_OPTION(msc_range_factor); GPO_LOAD_OPTION(msc_safety_factor); GPO_LOAD_OPTION(msc_lambda_limit); + GPO_LOAD_OPTION(msc_theta_limit); + GPO_LOAD_OPTION(angle_limit_factor); GPO_LOAD_OPTION(msc_step_algorithm); + GPO_LOAD_OPTION(form_factor); GPO_LOAD_OPTION(verbose); #undef GPO_LOAD_OPTION @@ -152,7 +168,10 @@ void to_json(nlohmann::json& j, GeantPhysicsOptions const& options) GPO_SAVE_OPTION(msc_range_factor); GPO_SAVE_OPTION(msc_safety_factor); GPO_SAVE_OPTION(msc_lambda_limit); + GPO_SAVE_OPTION(msc_theta_limit); + GPO_SAVE_OPTION(angle_limit_factor); GPO_SAVE_OPTION(msc_step_algorithm); + GPO_SAVE_OPTION(form_factor); GPO_SAVE_OPTION(verbose); #undef GPO_SAVE_OPTION diff --git a/src/celeritas/ext/detail/CelerEmStandardPhysics.cc b/src/celeritas/ext/detail/CelerEmStandardPhysics.cc index 9bc7eba8aa..5038ebd257 100644 --- a/src/celeritas/ext/detail/CelerEmStandardPhysics.cc +++ b/src/celeritas/ext/detail/CelerEmStandardPhysics.cc @@ -69,6 +69,28 @@ from_msc_step_algorithm(MscStepLimitAlgorithm const& msc_step_algorithm) } } +//---------------------------------------------------------------------------// +/*! + * Safely switch from NuclearFormFactorType to G4NuclearFormfactorType. + */ +G4NuclearFormfactorType +from_form_factor_type(NuclearFormFactorType const& form_factor) +{ + switch (form_factor) + { + case NuclearFormFactorType::none: + return G4NuclearFormfactorType::fNoneNF; + case NuclearFormFactorType::exponential: + return G4NuclearFormfactorType::fExponentialNF; + case NuclearFormFactorType::gaussian: + return G4NuclearFormfactorType::fGaussianNF; + case NuclearFormFactorType::flat: + return G4NuclearFormfactorType::fFlatNF; + default: + CELER_ASSERT_UNREACHABLE(); + } +} + //---------------------------------------------------------------------------// /*! * Construct with physics options. @@ -93,6 +115,8 @@ CelerEmStandardPhysics::CelerEmStandardPhysics(Options const& options) em_parameters.SetAuger(options.relaxation == RelaxationSelection::all); em_parameters.SetIntegral(options.integral_approach); em_parameters.SetLinearLossLimit(options.linear_loss_limit); + em_parameters.SetNuclearFormfactorType( + from_form_factor_type(options.form_factor)); em_parameters.SetMscStepLimitType( from_msc_step_algorithm(options.msc_step_algorithm)); em_parameters.SetMscRangeFactor(options.msc_range_factor); @@ -105,6 +129,7 @@ CelerEmStandardPhysics::CelerEmStandardPhysics(Options const& options) em_parameters.SetMscLambdaLimit( native_value_to(options.msc_lambda_limit).value()); #endif + em_parameters.SetMscThetaLimit(options.msc_theta_limit); em_parameters.SetLowestElectronEnergy( value_as(options.lowest_electron_energy) * CLHEP::MeV); diff --git a/src/celeritas/global/CoreParams.cc b/src/celeritas/global/CoreParams.cc index e8355b9687..1f6b9e6a40 100644 --- a/src/celeritas/global/CoreParams.cc +++ b/src/celeritas/global/CoreParams.cc @@ -23,6 +23,7 @@ #include "corecel/sys/MemRegistry.hh" #include "corecel/sys/ScopedMem.hh" #include "geocel/GeoParamsOutput.hh" +#include "celeritas/em/params/WentzelOKVIParams.hh" #include "celeritas/geo/GeoMaterialParams.hh" // IWYU pragma: keep #include "celeritas/geo/GeoParams.hh" // IWYU pragma: keep #include "celeritas/geo/detail/BoundaryAction.hh" @@ -88,6 +89,10 @@ build_params_refs(CoreParams::Input const& p, CoreScalars const& scalars) ref.rng = get_ref(*p.rng); ref.sim = get_ref(*p.sim); ref.init = get_ref(*p.init); + if (p.wentzel) + { + ref.wentzel = get_ref(*p.wentzel); + } CELER_ENSURE(ref); return ref; diff --git a/src/celeritas/global/CoreParams.hh b/src/celeritas/global/CoreParams.hh index ac1ebce6bc..2e39e9171a 100644 --- a/src/celeritas/global/CoreParams.hh +++ b/src/celeritas/global/CoreParams.hh @@ -23,9 +23,7 @@ namespace celeritas { //---------------------------------------------------------------------------// class ActionRegistry; -class AtomicRelaxationParams; class CutoffParams; -class FluctuationParams; class GeoMaterialParams; class MaterialParams; class OutputRegistry; @@ -33,6 +31,7 @@ class ParticleParams; class PhysicsParams; class SimParams; class TrackInitParams; +class WentzelOKVIParams; //---------------------------------------------------------------------------// /*! @@ -52,6 +51,7 @@ class CoreParams final : public ParamsDataInterface using SPConstRng = std::shared_ptr; using SPConstSim = std::shared_ptr; using SPConstTrackInit = std::shared_ptr; + using SPConstWentzelOKVI = std::shared_ptr; using SPActionRegistry = std::shared_ptr; using SPOutputRegistry = std::shared_ptr; @@ -72,6 +72,7 @@ class CoreParams final : public ParamsDataInterface SPConstRng rng; SPConstSim sim; SPConstTrackInit init; + SPConstWentzelOKVI wentzel; SPActionRegistry action_reg; SPOutputRegistry output_reg; @@ -106,6 +107,7 @@ class CoreParams final : public ParamsDataInterface SPConstRng const& rng() const { return input_.rng; } SPConstSim const& sim() const { return input_.sim; } SPConstTrackInit const& init() const { return input_.init; } + SPConstWentzelOKVI const& wentzel() const { return input_.wentzel; } SPActionRegistry const& action_reg() const { return input_.action_reg; } SPOutputRegistry const& output_reg() const { return input_.output_reg; } //!@} diff --git a/src/celeritas/global/CoreTrackData.hh b/src/celeritas/global/CoreTrackData.hh index 0fa21a26a2..46f2abac39 100644 --- a/src/celeritas/global/CoreTrackData.hh +++ b/src/celeritas/global/CoreTrackData.hh @@ -10,6 +10,7 @@ #include "corecel/Assert.hh" #include "corecel/data/Collection.hh" #include "celeritas/Types.hh" +#include "celeritas/em/data/WentzelOKVIData.hh" #include "celeritas/geo/GeoData.hh" #include "celeritas/geo/GeoMaterialData.hh" #include "celeritas/mat/MaterialData.hh" @@ -65,6 +66,7 @@ struct CoreParamsData RngParamsData rng; SimParamsData sim; TrackInitParamsData init; + WentzelOKVIData wentzel; CoreScalars scalars; @@ -89,6 +91,7 @@ struct CoreParamsData rng = other.rng; sim = other.sim; init = other.init; + wentzel = other.wentzel; scalars = other.scalars; return *this; } diff --git a/src/celeritas/io/ImportData.cc b/src/celeritas/io/ImportData.cc index 798c2565e8..10a8932e77 100644 --- a/src/celeritas/io/ImportData.cc +++ b/src/celeritas/io/ImportData.cc @@ -7,6 +7,8 @@ //---------------------------------------------------------------------------// #include "ImportData.hh" +#include + #include "corecel/Assert.hh" #include "corecel/io/Logger.hh" #include "celeritas/UnitTypes.hh" @@ -49,5 +51,36 @@ void convert_to_native(ImportData* data) CELER_ENSURE(data->units == units::NativeTraits::label()); } +//---------------------------------------------------------------------------// +/*! + * Whether an imported model of the given class is present. + */ +bool has_model(ImportData const& data, ImportModelClass model_class) +{ + for (ImportProcess const& process : data.processes) + { + for (ImportModel const& model : process.models) + { + if (model.model_class == model_class) + { + return true; + } + } + } + return false; +} + +//---------------------------------------------------------------------------// +/*! + * Whether an imported MSC model of the given class is present. + */ +bool has_msc_model(ImportData const& data, ImportModelClass model_class) +{ + return std::any_of( + data.msc_models.begin(), + data.msc_models.end(), + [&](ImportMscModel const& m) { return m.model_class == model_class; }); +} + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/io/ImportData.hh b/src/celeritas/io/ImportData.hh index dc3d3713ea..a56b53c1c9 100644 --- a/src/celeritas/io/ImportData.hh +++ b/src/celeritas/io/ImportData.hh @@ -88,5 +88,11 @@ struct ImportData // Recursively convert imported data to the native unit type void convert_to_native(ImportData* data); +// Whether an imported model of the given class is present +bool has_model(ImportData const&, ImportModelClass); + +// Whether an imported MSC model of the given class is present +bool has_msc_model(ImportData const&, ImportModelClass); + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/io/ImportParameters.hh b/src/celeritas/io/ImportParameters.hh index 7d03db176d..6dfdde7e2d 100644 --- a/src/celeritas/io/ImportParameters.hh +++ b/src/celeritas/io/ImportParameters.hh @@ -9,6 +9,7 @@ #include +#include "celeritas/Constants.hh" #include "celeritas/Types.hh" #include "celeritas/Units.hh" @@ -52,10 +53,16 @@ struct ImportEmParameters double msc_safety_factor{0.6}; //! MSC lambda limit [length] double msc_lambda_limit{1 * units::millimeter}; + //! Polar angle limii between single and multiple Coulomb scattering + double msc_theta_limit{constants::pi}; //! Kill secondaries below production cut bool apply_cuts{false}; //! Nuclear screening factor for single/multiple Coulomb scattering double screening_factor{1}; + //! Factor for dynamic computation of angular limit between SS and MSC + double angle_limit_factor{1}; + //! Nuclear form factor model for Coulomm scattering + NuclearFormFactorType form_factor{NuclearFormFactorType::exponential}; //! Whether parameters are assigned and valid explicit operator bool() const @@ -64,7 +71,9 @@ struct ImportEmParameters && msc_step_algorithm != MscStepLimitAlgorithm::size_ && msc_range_factor > 0 && msc_range_factor < 1 && msc_safety_factor >= 0.1 && msc_lambda_limit > 0 - && screening_factor > 0; + && msc_theta_limit >= 0 && msc_theta_limit <= constants::pi + && screening_factor > 0 && angle_limit_factor > 0 + && form_factor != NuclearFormFactorType::size_; } }; diff --git a/src/celeritas/phys/ProcessBuilder.cc b/src/celeritas/phys/ProcessBuilder.cc index 9beb3ee70d..2a787d0fc7 100644 --- a/src/celeritas/phys/ProcessBuilder.cc +++ b/src/celeritas/phys/ProcessBuilder.cc @@ -68,7 +68,6 @@ ProcessBuilder::ProcessBuilder(ImportData const& data, , brem_combined_(options.brem_combined) , enable_lpm_(data.em_params.lpm) , use_integral_xs_(data.em_params.integral_approach) - , coulomb_screening_factor_(data.em_params.screening_factor) { CELER_EXPECT(input_.material); CELER_EXPECT(input_.particle); @@ -241,12 +240,11 @@ auto ProcessBuilder::build_annihilation() -> SPProcess //---------------------------------------------------------------------------// auto ProcessBuilder::build_coulomb() -> SPProcess { - CoulombScatteringModel::Options options; - options.screening_factor = coulomb_screening_factor_; + CoulombScatteringProcess::Options options; options.use_integral_xs = use_integral_xs_; return std::make_shared( - this->particle(), this->material(), this->imported(), options); + this->particle(), this->imported(), options); } //---------------------------------------------------------------------------// diff --git a/src/celeritas/phys/ProcessBuilder.hh b/src/celeritas/phys/ProcessBuilder.hh index b2fed05681..392ed070b1 100644 --- a/src/celeritas/phys/ProcessBuilder.hh +++ b/src/celeritas/phys/ProcessBuilder.hh @@ -120,7 +120,6 @@ class ProcessBuilder bool brem_combined_; bool enable_lpm_; bool use_integral_xs_; - real_type coulomb_screening_factor_; //// HELPER FUNCTIONS //// diff --git a/test/celeritas/GlobalTestBase.cc b/test/celeritas/GlobalTestBase.cc index d5d113e26a..3e2949e974 100644 --- a/test/celeritas/GlobalTestBase.cc +++ b/test/celeritas/GlobalTestBase.cc @@ -94,6 +94,7 @@ auto GlobalTestBase::build_core() -> SPConstCore inp.rng = this->rng(); inp.sim = this->sim(); inp.init = this->init(); + inp.wentzel = this->wentzel(); inp.action_reg = this->action_reg(); inp.output_reg = this->output_reg(); CELER_ASSERT(inp); diff --git a/test/celeritas/GlobalTestBase.hh b/test/celeritas/GlobalTestBase.hh index e425f21226..b5a6639887 100644 --- a/test/celeritas/GlobalTestBase.hh +++ b/test/celeritas/GlobalTestBase.hh @@ -31,6 +31,7 @@ class ParticleParams; class PhysicsParams; class SimParams; class TrackInitParams; +class WentzelOKVIParams; class CoreParams; class OutputRegistry; @@ -67,6 +68,7 @@ class GlobalTestBase : public Test using SPConstRng = SP; using SPConstSim = SP; using SPConstTrackInit = SP; + using SPConstWentzelOKVI = SP; using SPConstCore = SP; using SPActionRegistry = SP; @@ -97,6 +99,7 @@ class GlobalTestBase : public Test inline SPConstRng const& rng(); inline SPConstSim const& sim(); inline SPConstTrackInit const& init(); + inline SPConstWentzelOKVI const& wentzel(); inline SPActionRegistry const& action_reg(); inline SPConstCore const& core(); inline SPConstCerenkov const& cerenkov(); @@ -113,6 +116,7 @@ class GlobalTestBase : public Test inline SPConstRng const& rng() const; inline SPConstSim const& sim() const; inline SPConstTrackInit const& init() const; + inline SPConstWentzelOKVI const& wentzel() const; inline SPActionRegistry const& action_reg() const; inline SPConstCore const& core() const; inline SPConstCerenkov const& cerenkov() const; @@ -138,6 +142,7 @@ class GlobalTestBase : public Test [[nodiscard]] virtual SPConstPhysics build_physics() = 0; [[nodiscard]] virtual SPConstSim build_sim() = 0; [[nodiscard]] virtual SPConstTrackInit build_init() = 0; + [[nodiscard]] virtual SPConstWentzelOKVI build_wentzel() = 0; [[nodiscard]] virtual SPConstAction build_along_step() = 0; [[nodiscard]] virtual SPConstCerenkov build_cerenkov() = 0; [[nodiscard]] virtual SPConstProperties build_properties() = 0; @@ -160,6 +165,7 @@ class GlobalTestBase : public Test SPConstRng rng_; SPConstSim sim_; SPConstTrackInit init_; + SPConstWentzelOKVI wentzel_; SPConstCore core_; SPOutputRegistry output_reg_; SPConstCerenkov cerenkov_; @@ -202,6 +208,18 @@ DEF_GTB_ACCESSORS(SPConstCore, core) DEF_GTB_ACCESSORS(SPConstCerenkov, cerenkov) DEF_GTB_ACCESSORS(SPConstProperties, properties) DEF_GTB_ACCESSORS(SPConstScintillation, scintillation) +auto GlobalTestBase::wentzel() -> SPConstWentzelOKVI const& +{ + if (!this->wentzel_) + { + this->wentzel_ = this->build_wentzel(); + } + return this->wentzel_; +} +auto GlobalTestBase::wentzel() const -> SPConstWentzelOKVI const& +{ + return this->wentzel_; +} #undef DEF_GTB_ACCESSORS diff --git a/test/celeritas/ImportedDataTestBase.cc b/test/celeritas/ImportedDataTestBase.cc index ae2e788cb7..c2f45c7e4b 100644 --- a/test/celeritas/ImportedDataTestBase.cc +++ b/test/celeritas/ImportedDataTestBase.cc @@ -7,6 +7,7 @@ //---------------------------------------------------------------------------// #include "ImportedDataTestBase.hh" +#include "celeritas/em/params/WentzelOKVIParams.hh" #include "celeritas/geo/GeoMaterialParams.hh" #include "celeritas/io/ImportData.hh" #include "celeritas/mat/MaterialParams.hh" @@ -70,6 +71,13 @@ auto ImportedDataTestBase::build_sim() -> SPConstSim return SimParams::from_import(this->imported_data(), this->particle()); } +//---------------------------------------------------------------------------// +auto ImportedDataTestBase::build_wentzel() -> SPConstWentzelOKVI +{ + return WentzelOKVIParams::from_import(this->imported_data(), + this->material()); +} + //---------------------------------------------------------------------------// auto ImportedDataTestBase::build_physics() -> SPConstPhysics { diff --git a/test/celeritas/ImportedDataTestBase.hh b/test/celeritas/ImportedDataTestBase.hh index a18c91cf44..fb10e6ae0d 100644 --- a/test/celeritas/ImportedDataTestBase.hh +++ b/test/celeritas/ImportedDataTestBase.hh @@ -50,6 +50,7 @@ class ImportedDataTestBase : virtual public GlobalGeoTestBase SPConstCutoff build_cutoff() override; SPConstPhysics build_physics() override; SPConstSim build_sim() override; + SPConstWentzelOKVI build_wentzel() override; SPConstCerenkov build_cerenkov() override; SPConstProperties build_properties() override; SPConstScintillation build_scintillation() override; diff --git a/test/celeritas/MockTestBase.hh b/test/celeritas/MockTestBase.hh index 9d82e7610c..7812a319b7 100644 --- a/test/celeritas/MockTestBase.hh +++ b/test/celeritas/MockTestBase.hh @@ -80,6 +80,7 @@ class MockTestBase : virtual public GlobalGeoTestBase, public OnlyCoreTestBase SPConstAction build_along_step() override; SPConstSim build_sim() override; SPConstTrackInit build_init() override; + SPConstWentzelOKVI build_wentzel() override { return nullptr; } virtual PhysicsOptions build_physics_options() const; diff --git a/test/celeritas/OnlyGeoTestBase.hh b/test/celeritas/OnlyGeoTestBase.hh index f8ed554a66..a17f44f798 100644 --- a/test/celeritas/OnlyGeoTestBase.hh +++ b/test/celeritas/OnlyGeoTestBase.hh @@ -33,6 +33,7 @@ class OnlyGeoTestBase : virtual public GlobalTestBase { CELER_ASSERT_UNREACHABLE(); } + SPConstWentzelOKVI build_wentzel() override { CELER_ASSERT_UNREACHABLE(); } }; //---------------------------------------------------------------------------// } // namespace test diff --git a/test/celeritas/SimpleTestBase.cc b/test/celeritas/SimpleTestBase.cc index 335e6c81ea..5f67db7075 100644 --- a/test/celeritas/SimpleTestBase.cc +++ b/test/celeritas/SimpleTestBase.cc @@ -8,6 +8,7 @@ #include "SimpleTestBase.hh" #include "celeritas/Quantities.hh" +#include "celeritas/em/params/WentzelOKVIParams.hh" #include "celeritas/em/process/ComptonProcess.hh" #include "celeritas/geo/GeoMaterialParams.hh" #include "celeritas/global/ActionRegistry.hh" @@ -181,6 +182,13 @@ auto SimpleTestBase::build_init() -> SPConstTrackInit return std::make_shared(input); } +//---------------------------------------------------------------------------// +auto SimpleTestBase::build_wentzel() -> SPConstWentzelOKVI +{ + WentzelOKVIParams::Options options; + return std::make_shared(this->material(), options); +} + //---------------------------------------------------------------------------// auto SimpleTestBase::build_along_step() -> SPConstAction { diff --git a/test/celeritas/SimpleTestBase.hh b/test/celeritas/SimpleTestBase.hh index 05589b1542..6fb81b3953 100644 --- a/test/celeritas/SimpleTestBase.hh +++ b/test/celeritas/SimpleTestBase.hh @@ -35,6 +35,7 @@ class SimpleTestBase : virtual public GlobalGeoTestBase, public OnlyCoreTestBase SPConstSim build_sim() override; SPConstTrackInit build_init() override; SPConstAction build_along_step() override; + SPConstWentzelOKVI build_wentzel() override; }; //---------------------------------------------------------------------------// diff --git a/test/celeritas/em/CoulombScattering.test.cc b/test/celeritas/em/CoulombScattering.test.cc index 6a75cf15be..bd28a14832 100644 --- a/test/celeritas/em/CoulombScattering.test.cc +++ b/test/celeritas/em/CoulombScattering.test.cc @@ -9,6 +9,7 @@ #include "celeritas/Units.hh" #include "celeritas/em/interactor/CoulombScatteringInteractor.hh" #include "celeritas/em/model/CoulombScatteringModel.hh" +#include "celeritas/em/params/WentzelOKVIParams.hh" #include "celeritas/em/process/CoulombScatteringProcess.hh" #include "celeritas/em/xs/WentzelTransportXsCalculator.hh" #include "celeritas/io/ImportParameters.hh" @@ -85,15 +86,15 @@ class CoulombScatteringTest : public InteractorHostTestBase {std::move(ip_electron), std::move(ip_positron)}); } - // Use default options - CoulombScatteringModel::Options options; + // Default to single scattering + WentzelOKVIParams::Options options; + options.is_combined = false; + options.polar_angle_limit = 0; + wentzel_ = std::make_shared(this->material_params(), + options); model_ = std::make_shared( - ActionId{0}, - *this->particle_params(), - *this->material_params(), - options, - this->imported_processes()); + ActionId{0}, *this->particle_params(), this->imported_processes()); // Set cutoffs CutoffParams::Input input; @@ -132,63 +133,113 @@ class CoulombScatteringTest : public InteractorHostTestBase } protected: + std::shared_ptr wentzel_; std::shared_ptr model_; + IsotopeComponentId isocomp_id_{0}; + ElementComponentId elcomp_id_{0}; + ElementId el_id_{0}; + MaterialId mat_id_{0}; }; -TEST_F(CoulombScatteringTest, wokvi_xs) +TEST_F(CoulombScatteringTest, helper) { - CoulombScatteringHostRef const& data = model_->host_ref(); - + struct Result + { + using VecReal = std::vector; + + VecReal screen_z; + VecReal scaled_kin_factor; + VecReal cos_thetamax_elec; + VecReal cos_thetamax_nuc; + VecReal xs_elec; + VecReal xs_nuc; + VecReal xs_ratio; + }; + + auto const material = this->material_track().make_material_view(); AtomicNumber const target_z - = this->material_params()->get(ElementId{0}).atomic_number(); + = this->material_params()->get(el_id_).atomic_number(); - MevEnergy const cutoff = this->cutoff_params() - ->get(MaterialId{0}) - .energy(this->particle_track().particle_id()); + MevEnergy const cutoff = this->cutoff_params()->get(mat_id_).energy( + this->particle_track().particle_id()); - std::vector const energies = {50, 100, 200, 1000, 13000}; + real_type const cos_thetamax = model_->host_ref().cos_thetamax(); - static real_type const expected_screen_z[] = {2.1181757502465e-08, - 5.3641196710457e-09, - 1.3498490873627e-09, - 5.4280909096648e-11, - 3.2158426877075e-13}; + Result result; + for (real_type energy : {50, 100, 200, 1000, 13000}) + { + this->set_inc_particle(pdg::electron(), MevEnergy{energy}); - static real_type const expected_cos_t_max[] = {0.99989885103277, - 0.99997458240728, - 0.99999362912075, - 0.99999974463379, - 0.99999999848823}; + WentzelHelper helper(this->particle_track(), + material, + target_z, + wentzel_->host_ref(), + model_->host_ref().ids, + cutoff); + + EXPECT_SOFT_EQ(1.1682, helper.mott_factor()); + result.screen_z.push_back(helper.screening_coefficient()); + // Scale the xs factor by 1 / r_e^2 so the values will be large enough + // for the soft equivalence comparison to catch any differences + result.scaled_kin_factor.push_back(helper.kin_factor() + / ipow<2>(constants::r_electron)); + result.cos_thetamax_elec.push_back(helper.cos_thetamax_electron()); + real_type const cos_thetamax_nuc = helper.cos_thetamax_nuclear(); + result.cos_thetamax_nuc.push_back(cos_thetamax_nuc); + result.xs_elec.push_back( + helper.calc_xs_electron(cos_thetamax_nuc, cos_thetamax) + / units::barn); + result.xs_nuc.push_back( + helper.calc_xs_nuclear(cos_thetamax_nuc, cos_thetamax) + / units::barn); + result.xs_ratio.push_back( + helper.calc_xs_ratio(cos_thetamax_nuc, cos_thetamax)); + } - static real_type const expected_xsecs[] = {0.033319844069031, + static double const expected_screen_z[] = {2.1181757502465e-08, + 5.3641196710457e-09, + 1.3498490873627e-09, + 5.4280909096648e-11, + 3.2158426877075e-13}; + static double const expected_scaled_kin_factor[] = {0.018652406309778, + 0.0047099159161888, + 0.0011834423911797, + 4.7530717872407e-05, + 2.8151208086621e-07}; + static double const expected_cos_thetamax_elec[] = {0.99989885103277, + 0.99997458240728, + 0.99999362912075, + 0.99999974463379, + 0.99999999848823}; + static double const expected_cos_thetamax_nuc[] = {1, 1, 1, 1, 1}; + static double const expected_xs_elec[] = {40826.46816866, + 40708.229862005, + 40647.018860182, + 40596.955206725, + 40585.257368735}; + static double const expected_xs_nuc[] = {1184463.4246675, + 1181036.9405781, + 1179263.0534548, + 1177812.2021574, + 1177473.1955599}; + static double const expected_xs_ratio[] = {0.033319844069031, 0.033319738720425, 0.033319684608429, 0.033319640583261, 0.03331963032739}; - - std::vector xsecs, cos_t_maxs, screen_zs; - for (real_type energy : energies) - { - this->set_inc_particle(pdg::electron(), MevEnergy{energy}); - - WentzelHelper helper(this->particle_track(), target_z, data, cutoff); - - xsecs.push_back(helper.calc_xs_ratio()); - cos_t_maxs.push_back(helper.costheta_max_electron()); - screen_zs.push_back(helper.screening_coefficient()); - } - - EXPECT_VEC_SOFT_EQ(expected_xsecs, xsecs); - EXPECT_VEC_SOFT_EQ(expected_screen_z, screen_zs); - EXPECT_VEC_SOFT_EQ(expected_cos_t_max, cos_t_maxs); + EXPECT_VEC_SOFT_EQ(expected_screen_z, result.screen_z); + EXPECT_VEC_SOFT_EQ(expected_scaled_kin_factor, result.scaled_kin_factor); + EXPECT_VEC_SOFT_EQ(expected_cos_thetamax_elec, result.cos_thetamax_elec); + EXPECT_VEC_SOFT_EQ(expected_cos_thetamax_nuc, result.cos_thetamax_nuc); + EXPECT_VEC_SOFT_EQ(expected_xs_elec, result.xs_elec); + EXPECT_VEC_SOFT_EQ(expected_xs_nuc, result.xs_nuc); + EXPECT_VEC_SOFT_EQ(expected_xs_ratio, result.xs_ratio); } TEST_F(CoulombScatteringTest, mott_xs) { - CoulombScatteringHostRef const& data = model_->host_ref(); - - CoulombScatteringElementData const& element_data - = data.elem_data[ElementId(0)]; + MottElementData const& element_data + = wentzel_->host_ref().elem_data[el_id_]; MottRatioCalculator xsec(element_data, sqrt(this->particle_track().beta_sq())); @@ -216,14 +267,11 @@ TEST_F(CoulombScatteringTest, mott_xs) TEST_F(CoulombScatteringTest, wokvi_transport_xs) { - // Copper - MaterialId mat_id{0}; - ElementId elm_id{0}; - - AtomicNumber const z = this->material_params()->get(elm_id).atomic_number(); + auto const material = this->material_track().make_material_view(); + AtomicNumber const z = this->material_params()->get(el_id_).atomic_number(); // Incident particle energy cutoff - MevEnergy const cutoff = this->cutoff_params()->get(mat_id).energy( + MevEnergy const cutoff = this->cutoff_params()->get(mat_id_).energy( this->particle_track().particle_id()); std::vector xs; @@ -232,13 +280,18 @@ TEST_F(CoulombScatteringTest, wokvi_transport_xs) this->set_inc_particle(pdg::electron(), MevEnergy{energy}); auto const& particle = this->particle_track(); - WentzelHelper helper(particle, z, model_->host_ref(), cutoff); + WentzelHelper helper(particle, + material, + z, + wentzel_->host_ref(), + model_->host_ref().ids, + cutoff); WentzelTransportXsCalculator calc_transport_xs(particle, helper); - for (real_type costheta_max : {-1.0, -0.5, 0.0, 0.5, 0.75, 0.99, 1.0}) + for (real_type cos_thetamax : {-1.0, -0.5, 0.0, 0.5, 0.75, 0.99, 1.0}) { // Get cross section in barns - xs.push_back(calc_transport_xs(costheta_max) / units::barn); + xs.push_back(calc_transport_xs(cos_thetamax) / units::barn); } } static double const expected_xs[] = {0.18738907324438, @@ -283,11 +336,10 @@ TEST_F(CoulombScatteringTest, simple_scattering) { int const num_samples = 4; - IsotopeView const isotope = this->material_track() - .make_material_view() - .make_element_view(ElementComponentId{0}) - .make_isotope_view(IsotopeComponentId{0}); - auto cutoffs = this->cutoff_params()->get(MaterialId{0}); + auto const material = this->material_track().make_material_view(); + IsotopeView const isotope + = material.make_element_view(elcomp_id_).make_isotope_view(isocomp_id_); + auto cutoffs = this->cutoff_params()->get(mat_id_); RandomEngine& rng_engine = this->rng(); @@ -299,10 +351,12 @@ TEST_F(CoulombScatteringTest, simple_scattering) { this->set_inc_particle(pdg::electron(), MevEnergy{energy}); CoulombScatteringInteractor interact(model_->host_ref(), + wentzel_->host_ref(), this->particle_track(), this->direction(), + material, isotope, - ElementId{0}, + el_id_, cutoffs); for ([[maybe_unused]] int i : range(num_samples)) @@ -317,52 +371,52 @@ TEST_F(CoulombScatteringTest, simple_scattering) } } static double const expected_cos_theta[] = {1, - 0.9996601312603, - 0.99998628518524, - 0.9998960794451, - 1, - 0.99991152273528, - 0.99998247385882, - 0.99988427878015, + 0.99950360343422, + 0.98776892641281, + 0.99837727448607, 1, - 0.99999970191006, - 0.99999637934855, - 0.99999885954183, - 0.999999997443, - 0.99999999406847, - 0.99999999849197, + 0.9999716884097, + 0.99985707764428, + 0.99997835395879, 1, - 0.99999999992169, + 0.99999688465904, + 0.99999974351257, + 0.99999918571981, + 0.99999995498814, + 0.99999998059604, + 0.99999992367847, 1, - 0.99999999209019, - 0.99999999999489, + 0.99999999984949, 1, - 0.99999999999999, + 0.99999999999851, + 0.99999999769513, + 0.99999999999996, + 0.99999999999998, 0.99999999999999, 1}; static double const expected_delta_energy[] = {0, - 1.4170232209842e-09, - 5.7181509527382e-11, - 4.3327857968123e-10, + 2.069638599389e-09, + 5.0995313499724e-08, + 6.7656699409557e-09, 0, - 3.0519519134131e-09, - 6.0455007666604e-10, - 3.9917101846143e-09, + 9.7658547915103e-10, + 4.9299914151035e-09, + 7.4666273164326e-10, 0, - 5.6049742624964e-10, - 6.8078875870015e-09, - 2.1443966602419e-09, - 4.406643938637e-10, - 1.0222294122286e-09, - 2.5988811103161e-10, + 5.8577551698136e-09, + 4.8227200011297e-10, + 1.5310863688001e-09, + 7.7572508416779e-09, + 3.3440414881625e-09, + 1.315311237704e-08, 0, - 1.3371845852816e-09, + 2.5702320272103e-09, 0, - 1.350750835627e-07, - 8.7311491370201e-11, - 4.8021320253611e-10, - 8.7311491370201e-10, - 1.6734702512622e-09, + 2.5465851649642e-11, + 3.9359974834952e-08, + 7.1158865466714e-09, + 3.9435690268874e-09, + 2.0663719624281e-09, 0}; EXPECT_VEC_SOFT_EQ(expected_cos_theta, cos_theta); EXPECT_VEC_SOFT_EQ(expected_delta_energy, delta_energy); @@ -370,28 +424,33 @@ TEST_F(CoulombScatteringTest, simple_scattering) TEST_F(CoulombScatteringTest, distribution) { - CoulombScatteringHostRef const& data = model_->host_ref(); + auto const material = this->material_track().make_material_view(); + IsotopeView const isotope + = material.make_element_view(elcomp_id_).make_isotope_view(isocomp_id_); + + // TODO: Use proton ParticleId{2} + MevEnergy const cutoff + = this->cutoff_params()->get(mat_id_).energy(ParticleId{0}); + std::vector avg_angles; for (real_type energy : {1, 50, 100, 200, 1000, 13000}) { this->set_inc_particle(pdg::electron(), MevEnergy{energy}); - CoulombScatteringElementData const& element_data - = data.elem_data[ElementId(0)]; - - IsotopeView const isotope - = this->material_track() - .make_material_view() - .make_element_view(ElementComponentId{0}) - .make_isotope_view(IsotopeComponentId{0}); - - // TODO: Use proton ParticleId{2} - MevEnergy const cutoff - = this->cutoff_params()->get(MaterialId{0}).energy(ParticleId{0}); - - WentzelDistribution distrib( - this->particle_track(), isotope, element_data, cutoff, data); + WentzelHelper helper(this->particle_track(), + material, + isotope.atomic_number(), + wentzel_->host_ref(), + model_->host_ref().ids, + cutoff); + WentzelDistribution sample_angle(wentzel_->host_ref(), + helper, + this->particle_track(), + isotope, + el_id_, + helper.cos_thetamax_nuclear(), + model_->host_ref().cos_thetamax()); RandomEngine& rng_engine = this->rng(); @@ -400,19 +459,19 @@ TEST_F(CoulombScatteringTest, distribution) int const num_samples = 4096; for ([[maybe_unused]] int i : range(num_samples)) { - avg_angle += distrib(rng_engine); + avg_angle += sample_angle(rng_engine); } avg_angle /= num_samples; avg_angles.push_back(avg_angle); } - static double const expected_avg_angles[] = {0.99933909299229, - 0.99999960697043, - 0.99999986035881, - 0.99999998052954, - 0.99999999917037, - 0.9999999999969}; + static double const expected_avg_angles[] = {0.99957853627426, + 0.99999954645904, + 0.99999989882947, + 0.99999996985799, + 0.99999999945722, + 0.99999999999487}; EXPECT_VEC_SOFT_EQ(expected_avg_angles, avg_angles); } diff --git a/test/celeritas/ext/GeantImporter.test.cc b/test/celeritas/ext/GeantImporter.test.cc index d14b6f4ff7..773862826d 100644 --- a/test/celeritas/ext/GeantImporter.test.cc +++ b/test/celeritas/ext/GeantImporter.test.cc @@ -256,7 +256,7 @@ class FourSteelSlabsEmStandard : public GeantImporterTest { nlohmann::json out = opts; EXPECT_JSON_EQ( - R"json({"annihilation":true,"apply_cuts":false,"brems":"all","compton_scattering":true,"coulomb_scattering":false,"default_cutoff":0.1,"eloss_fluctuation":true,"em_bins_per_decade":7,"gamma_conversion":true,"gamma_general":false,"integral_approach":true,"ionization":true,"linear_loss_limit":0.01,"lowest_electron_energy":[0.001,"MeV"],"lpm":true,"max_energy":[100000000.0,"MeV"],"min_energy":[0.0001,"MeV"],"msc":"urban","msc_lambda_limit":0.1,"msc_range_factor":0.04,"msc_safety_factor":0.6,"msc_step_algorithm":"safety","photoelectric":true,"rayleigh_scattering":true,"relaxation":"all","verbose":true})json", + R"json({"angle_limit_factor":1.0,"annihilation":true,"apply_cuts":false,"brems":"all","compton_scattering":true,"coulomb_scattering":false,"default_cutoff":0.1,"eloss_fluctuation":true,"em_bins_per_decade":7,"form_factor":"exponential","gamma_conversion":true,"gamma_general":false,"integral_approach":true,"ionization":true,"linear_loss_limit":0.01,"lowest_electron_energy":[0.001,"MeV"],"lpm":true,"max_energy":[100000000.0,"MeV"],"min_energy":[0.0001,"MeV"],"msc":"urban","msc_lambda_limit":0.1,"msc_range_factor":0.04,"msc_safety_factor":0.6,"msc_step_algorithm":"safety","msc_theta_limit":3.141592653589793,"photoelectric":true,"rayleigh_scattering":true,"relaxation":"all","verbose":true})json", std::string(out.dump())); } #endif @@ -805,9 +805,14 @@ TEST_F(FourSteelSlabsEmStandard, em_parameters) EXPECT_DOUBLE_EQ(0.01, em_params.linear_loss_limit); EXPECT_DOUBLE_EQ(0.001, em_params.lowest_electron_energy); EXPECT_EQ(true, em_params.auger); + EXPECT_EQ(MscStepLimitAlgorithm::safety, em_params.msc_step_algorithm); EXPECT_DOUBLE_EQ(0.04, em_params.msc_range_factor); EXPECT_DOUBLE_EQ(0.6, em_params.msc_safety_factor); EXPECT_REAL_EQ(0.1, to_cm(em_params.msc_lambda_limit)); + EXPECT_REAL_EQ(constants::pi, em_params.msc_theta_limit); + EXPECT_EQ(false, em_params.apply_cuts); + EXPECT_EQ(1, em_params.screening_factor); + EXPECT_EQ(1, em_params.angle_limit_factor); } //---------------------------------------------------------------------------// diff --git a/test/celeritas/mat/Material.test.cc b/test/celeritas/mat/Material.test.cc index 290fd34b89..cbd72882b5 100644 --- a/test/celeritas/mat/Material.test.cc +++ b/test/celeritas/mat/Material.test.cc @@ -195,17 +195,17 @@ TEST_F(MaterialTest, material_view) { // vacuum MaterialView mat = params->get(MaterialId{1}); - EXPECT_SOFT_EQ(0, mat.number_density()); - EXPECT_SOFT_EQ(0, mat.temperature()); + EXPECT_EQ(0, mat.number_density()); + EXPECT_EQ(0, mat.temperature()); EXPECT_EQ(MatterState::unspecified, mat.matter_state()); - EXPECT_SOFT_EQ(0, mat.zeff()); - EXPECT_SOFT_EQ(0, mat.density()); - EXPECT_SOFT_EQ(0, mat.electron_density()); - EXPECT_SOFT_EQ(std::numeric_limits::infinity(), - mat.radiation_length()); - EXPECT_SOFT_EQ(0, mat.mean_excitation_energy().value()); - EXPECT_SOFT_EQ(-std::numeric_limits::infinity(), - mat.log_mean_excitation_energy().value()); + EXPECT_EQ(0, mat.zeff()); + EXPECT_EQ(0, mat.density()); + EXPECT_EQ(0, mat.electron_density()); + EXPECT_EQ(std::numeric_limits::infinity(), + mat.radiation_length()); + EXPECT_EQ(0, mat.mean_excitation_energy().value()); + EXPECT_EQ(-std::numeric_limits::infinity(), + mat.log_mean_excitation_energy().value()); // Test element view auto els = mat.elements();