diff --git a/src/celeritas/neutron/data/NeutronElasticData.hh b/src/celeritas/neutron/data/NeutronElasticData.hh index 328a555abf..f65ca5aba4 100644 --- a/src/celeritas/neutron/data/NeutronElasticData.hh +++ b/src/celeritas/neutron/data/NeutronElasticData.hh @@ -147,7 +147,7 @@ struct NeutronElasticData NeutronElasticXsData xs; //! Neutron elastic differential cross section coefficients - NeutronElasticDiffXsData Q2; + NeutronElasticDiffXsData diff_xs; //// MEMBER FUNCTIONS //// @@ -176,7 +176,7 @@ struct NeutronElasticData ids = other.ids; neutron_mass = other.neutron_mass; xs = other.xs; - Q2 = other.Q2; + diff_xs = other.diff_xs; return *this; } }; diff --git a/src/celeritas/neutron/interactor/NeutronElasticInteractor.hh b/src/celeritas/neutron/interactor/NeutronElasticInteractor.hh index 23fecf0b9a..593e3e2d04 100644 --- a/src/celeritas/neutron/interactor/NeutronElasticInteractor.hh +++ b/src/celeritas/neutron/interactor/NeutronElasticInteractor.hh @@ -38,6 +38,7 @@ class NeutronElasticInteractor //! \name Type aliases using Energy = units::MevEnergy; using Mass = units::MevMass; + using Momentum = units::MevMomentum; //!@} public: @@ -69,8 +70,8 @@ class NeutronElasticInteractor IsotopeView const& target_; // Value of neutron MeVMass - real_type neutron_mass_; - real_type neutron_p_; + Mass neutron_mass_; + Momentum neutron_p_; // Sampler UniformRealDist sample_phi_; @@ -92,16 +93,14 @@ CELER_FUNCTION NeutronElasticInteractor::NeutronElasticInteractor( , inc_energy_(particle.energy()) , inc_direction_(inc_direction) , target_(target) - , neutron_mass_(value_as(shared_.neutron_mass)) - , neutron_p_(std::sqrt(inc_energy_.value() - * (inc_energy_.value() + 2 * neutron_mass_))) + , neutron_mass_(shared_.neutron_mass) + , neutron_p_(particle.momentum()) , sample_phi_(0, 2 * constants::pi) , sample_Q2_(shared_, target_, neutron_p_) { CELER_EXPECT(particle.particle_id() == shared_.ids.neutron); CELER_EXPECT(inc_energy_ > zero_quantity()); } - //---------------------------------------------------------------------------// /*! * Sample the final state of the neutron-nucleus elastic scattering. @@ -109,16 +108,16 @@ CELER_FUNCTION NeutronElasticInteractor::NeutronElasticInteractor( template CELER_FUNCTION Interaction NeutronElasticInteractor::operator()(Engine& rng) { - using Energy = units::MevEnergy; - // Scattered neutron with respect to the axis of incident direction Interaction result; // The momentum magnitude in the c.m. frame real_type target_mass = value_as(target_.nuclear_mass()); - real_type neutron_energy = neutron_mass_ + inc_energy_.value(); - real_type cm_p = neutron_p_ - / std::sqrt(1 + ipow<2>(neutron_mass_ / target_mass) + real_type neutron_mass = value_as(neutron_mass_); + real_type neutron_energy = neutron_mass + inc_energy_.value(); + + real_type cm_p = value_as(neutron_p_) + / std::sqrt(1 + ipow<2>(neutron_mass / target_mass) + 2 * neutron_energy / target_mass); // Sample the scattered direction from the invariant momentum transfer @@ -128,16 +127,18 @@ CELER_FUNCTION Interaction NeutronElasticInteractor::operator()(Engine& rng) // Boost to the center of mass (c.m.) frame Real3 cm_mom = cm_p * from_spherical(cos_theta, sample_phi_(rng)); - FourVector nlv1( - {cm_mom, std::sqrt(ipow<2>(cm_p) + ipow<2>(neutron_mass_))}); + FourVector nlv1({cm_mom, std::sqrt(ipow<2>(cm_p) + ipow<2>(neutron_mass))}); - FourVector lv({{0, 0, neutron_p_}, neutron_energy + target_mass}); + FourVector lv({{0, 0, value_as(neutron_p_)}, + neutron_energy + target_mass}); boost(boost_vector(lv), &nlv1); - result.direction = make_unit_vector(rotate(nlv1.mom, inc_direction_)); + + result.direction + = make_unit_vector(rotate(make_unit_vector(nlv1.mom), inc_direction_)); // Kinetic energy of the scattered neutron and the recoiled nucleus lv.energy -= nlv1.energy; - result.energy = Energy(nlv1.energy - neutron_mass_); + result.energy = Energy(nlv1.energy - neutron_mass); real_type recoil_energy = clamp_to_nonneg(lv.energy - target_mass); result.energy_deposition = Energy{recoil_energy}; diff --git a/src/celeritas/neutron/interactor/detail/InvariantQ2Sampler.hh b/src/celeritas/neutron/interactor/detail/InvariantQ2Sampler.hh index 1d04b7129d..80c665860e 100644 --- a/src/celeritas/neutron/interactor/detail/InvariantQ2Sampler.hh +++ b/src/celeritas/neutron/interactor/detail/InvariantQ2Sampler.hh @@ -38,13 +38,15 @@ class InvariantQ2Sampler //!@{ //! \name Type aliases using Mass = units::MevMass; + using Momentum = units::MevMomentum; + using AtomicMassNumber = AtomicNumber; //!@} public: // Construct with shared and target data, and the neutron momentum inline CELER_FUNCTION InvariantQ2Sampler(NeutronElasticRef const& shared, IsotopeView const& target, - real_type neutron_p); + Momentum neutron_p); // Sample the momentum transfer template @@ -59,19 +61,19 @@ class InvariantQ2Sampler ChipsDiffXsCoefficients::ChipsArray const& par_; - // Mass values in GeV - real_type neutron_mass_; - real_type target_mass_; + // Mass of neutron and target + Mass neutron_mass_; + Mass target_mass_; // Atomic mass number (A) of the target - real_type A_; + AtomicMassNumber A_; bool heavy_target_{false}; // Momentum magnitude of the incident neutron (GeV/c) - real_type neutron_p_; + Momentum neutron_p_; // Maximum momentum transfer for the elastic scattering - real_type max_Q2_; + real_type max_q2_; // Parameters for the CHIPS differential cross section - ExchangeParameters par_Q2_; + ExchangeParameters par_q2_; //// HELPER FUNCTIONS //// @@ -83,10 +85,10 @@ class InvariantQ2Sampler } // Calculate the maximum momentum transfer - inline CELER_FUNCTION real_type calc_max_Q2(real_type p) const; + inline CELER_FUNCTION real_type calc_max_q2(Momentum p) const; // Calculate parameters used in the quark-exchange process - inline CELER_FUNCTION ExchangeParameters calc_par_Q2(real_type p) const; + inline CELER_FUNCTION ExchangeParameters calc_par_q2(Momentum p) const; //// COMMON PROPERTIES //// @@ -94,7 +96,10 @@ class InvariantQ2Sampler static CELER_CONSTEXPR_FUNCTION real_type to_gev() { return 1e-3; } // S-wave for neutron p < 14 MeV/c (kinetic energy < 0.1 MeV) - static CELER_CONSTEXPR_FUNCTION real_type s_wave_limit() { return 1.4e-2; } + static CELER_CONSTEXPR_FUNCTION Momentum s_wave_limit() + { + return Momentum{14}; + } }; //---------------------------------------------------------------------------// @@ -109,17 +114,16 @@ class InvariantQ2Sampler CELER_FUNCTION InvariantQ2Sampler::InvariantQ2Sampler(NeutronElasticRef const& shared, IsotopeView const& target, - real_type neutron_p) - - : par_(shared.Q2.coeffs[target.isotope_id()].par) - , neutron_mass_(value_as(shared.neutron_mass) * this->to_gev()) - , target_mass_(value_as(target.nuclear_mass()) * this->to_gev()) - , A_(static_cast(target.atomic_mass_number().get())) - , heavy_target_(A_ > 6) - , neutron_p_(neutron_p * this->to_gev()) + Momentum neutron_p) + : par_(shared.diff_xs.coeffs[target.isotope_id()].par) + , neutron_mass_(shared.neutron_mass) + , target_mass_(target.nuclear_mass()) + , A_(target.atomic_mass_number()) + , heavy_target_(A_.get() > 6) + , neutron_p_(neutron_p) { - max_Q2_ = calc_max_Q2(neutron_p_); - par_Q2_ = calc_par_Q2(neutron_p_); + max_q2_ = calc_max_q2(neutron_p_); + par_q2_ = calc_par_q2(neutron_p_); } //---------------------------------------------------------------------------// @@ -135,25 +139,25 @@ CELER_FUNCTION auto InvariantQ2Sampler::operator()(Engine& rng) -> real_type // Sample \f$ Q^{2} \f$ below S-wave limit if (neutron_p_ < this->s_wave_limit()) { - return max_Q2_ * generate_canonical(rng) / ipow<2>(this->to_gev()); + return max_q2_ * generate_canonical(rng) / ipow<2>(this->to_gev()); } // Sample \f$ Q^{2} \f$ real_type q2 = 0; - if (A_ == 1) + if (A_ == AtomicMassNumber{1}) { // Special case for the \f$ n + p \rightarrow n + p \f$ channel - real_type R1 = -std::expm1(-max_Q2_ * par_Q2_.slope[0]); - real_type R2 = -std::expm1(-max_Q2_ * par_Q2_.slope[1]); + real_type R1 = -std::expm1(-max_q2_ * par_q2_.slope[0]); + real_type R2 = -std::expm1(-max_q2_ * par_q2_.slope[1]); - real_type I1 = R1 * par_Q2_.expnt[0]; - real_type I2 = R2 * par_Q2_.expnt[1] / par_Q2_.slope[1]; + real_type I1 = R1 * par_q2_.expnt[0]; + real_type I2 = R2 * par_q2_.expnt[1] / par_q2_.slope[1]; // Sample by t-channel and u-channel (charge exchange) q2 = (BernoulliDistribution(I1 / (I1 + I2))(rng)) - ? sample_q2(R1, rng) / par_Q2_.slope[0] - : max_Q2_ - sample_q2(R2, rng) / par_Q2_.slope[1]; + ? sample_q2(R1, rng) / par_q2_.slope[0] + : max_q2_ - sample_q2(R2, rng) / par_q2_.slope[1]; } else { @@ -163,46 +167,46 @@ CELER_FUNCTION auto InvariantQ2Sampler::operator()(Engine& rng) -> real_type constexpr real_type one_seventh = 1 / real_type(7); real_type R1 = -std::expm1( - -max_Q2_ * (par_Q2_.slope[0] + max_Q2_ * par_Q2_.ss)); + -max_q2_ * (par_q2_.slope[0] + max_q2_ * par_q2_.ss)); // power 3 for low-A and power 5 for high-A - real_type factor = heavy_target_ ? ipow<5>(max_Q2_) : ipow<3>(max_Q2_); - real_type R2 = -std::expm1(-factor * par_Q2_.slope[1]); + real_type factor = heavy_target_ ? ipow<5>(max_q2_) : ipow<3>(max_q2_); + real_type R2 = -std::expm1(-factor * par_q2_.slope[1]); // power 1 for low-A and power 7 for high-A - factor = heavy_target_ ? ipow<7>(max_Q2_) : max_Q2_; - real_type R3 = -std::expm1(-factor * par_Q2_.slope[2]); - real_type R4 = -std::expm1(-max_Q2_ * par_Q2_.slope[3]); - - real_type I1 = R1 * par_Q2_.expnt[0]; - real_type I2 = R2 * par_Q2_.expnt[1]; - real_type I3 = R3 * par_Q2_.expnt[2]; - real_type I4 = R4 * par_Q2_.expnt[3]; + factor = heavy_target_ ? ipow<7>(max_q2_) : max_q2_; + real_type R3 = -std::expm1(-factor * par_q2_.slope[2]); + real_type R4 = -std::expm1(-max_q2_ * par_q2_.slope[3]); + + real_type I1 = R1 * par_q2_.expnt[0]; + real_type I2 = R2 * par_q2_.expnt[1]; + real_type I3 = R3 * par_q2_.expnt[2]; + real_type I4 = R4 * par_q2_.expnt[3]; real_type I5 = I1 + I2; real_type I6 = I3 + I5; real_type rand = (I4 + I6) * generate_canonical(rng); if (rand < I1) { - real_type tss = 2 * par_Q2_.ss; - q2 = this->sample_q2(R1, rng) / par_Q2_.slope[0]; + real_type tss = 2 * par_q2_.ss; + q2 = this->sample_q2(R1, rng) / par_q2_.slope[0]; if (std::fabs(tss) > 1e-7) { - q2 = (std::sqrt(par_Q2_.slope[0] - * (par_Q2_.slope[0] + 2 * tss * q2)) - - par_Q2_.slope[0]) + q2 = (std::sqrt(par_q2_.slope[0] + * (par_q2_.slope[0] + 2 * tss * q2)) + - par_q2_.slope[0]) / tss; } } else if (rand < I5) { - q2 = clamp_to_nonneg(this->sample_q2(R2, rng) / par_Q2_.slope[1]); + q2 = clamp_to_nonneg(this->sample_q2(R2, rng) / par_q2_.slope[1]); q2 = (heavy_target_) ? std::pow(q2, one_fifth) : std::pow(q2, one_third); } else if (rand < I6) { - q2 = clamp_to_nonneg(this->sample_q2(R3, rng) / par_Q2_.slope[2]); + q2 = clamp_to_nonneg(this->sample_q2(R3, rng) / par_q2_.slope[2]); if (heavy_target_) { q2 = std::pow(q2, one_seventh); @@ -210,15 +214,15 @@ CELER_FUNCTION auto InvariantQ2Sampler::operator()(Engine& rng) -> real_type } else { - q2 = this->sample_q2(R4, rng) / par_Q2_.slope[3]; + q2 = this->sample_q2(R4, rng) / par_q2_.slope[3]; // u reduced for light-A (starts from 0) if (!heavy_target_) { - q2 = max_Q2_ - q2; + q2 = max_q2_ - q2; } } } - return clamp(q2, real_type{0}, max_Q2_) / ipow<2>(this->to_gev()); + return clamp(q2, real_type{0}, max_q2_) / ipow<2>(this->to_gev()); } //---------------------------------------------------------------------------// @@ -235,15 +239,18 @@ CELER_FUNCTION auto InvariantQ2Sampler::operator()(Engine& rng) -> real_type * collision angle is 90 degree in the center of mass system, is currently * excluded, but may be supported if there is a user case. */ -CELER_FUNCTION real_type InvariantQ2Sampler::calc_max_Q2(real_type p) const +CELER_FUNCTION real_type InvariantQ2Sampler::calc_max_q2(Momentum p) const { // Momentum and mass square of the incident neutron - real_type p2 = ipow<2>(p); - real_type m2 = ipow<2>(neutron_mass_); - real_type a2 = ipow<2>(target_mass_); + real_type target_mass = value_as(target_mass_); + + real_type p2 = ipow<2>(value_as(p)); + real_type m2 = ipow<2>(value_as(neutron_mass_)); + real_type a2 = ipow<2>(target_mass); // Return the maximum momentum transfer in the value of clhep::GeV square - return 4 * a2 * p2 / (2 * target_mass_ * std::sqrt(p2 + m2) + m2 + a2); + return 4 * a2 * p2 * ipow<2>(this->to_gev()) + / (2 * target_mass * std::sqrt(p2 + m2) + m2 + a2); } //---------------------------------------------------------------------------// @@ -253,8 +260,11 @@ CELER_FUNCTION real_type InvariantQ2Sampler::calc_max_Q2(real_type p) const * \param p the neutron momentum in the lab frame (value in clhep::GeV unit). */ CELER_FUNCTION -auto InvariantQ2Sampler::calc_par_Q2(real_type p) const -> ExchangeParameters +auto InvariantQ2Sampler::calc_par_q2(Momentum neutron_p) const + -> ExchangeParameters { + // ExchangeParameters + real_type p = value_as(neutron_p) * this->to_gev(); real_type lp = std::log(p); real_type sp = std::sqrt(p); real_type p2 = ipow<2>(p); @@ -263,7 +273,7 @@ auto InvariantQ2Sampler::calc_par_Q2(real_type p) const -> ExchangeParameters ExchangeParameters result; - if (A_ == 1) + if (A_ == AtomicMassNumber{1}) { // Special case for the \f$ n + p \rightarrow n + p \f$ channel real_type dl1 = lp - par_[3]; @@ -305,7 +315,7 @@ auto InvariantQ2Sampler::calc_par_Q2(real_type p) const -> ExchangeParameters if (!heavy_target_) { - real_type pah = std::pow(p, real_type(0.5) * A_); + real_type pah = std::pow(p, real_type(0.5) * A_.get()); real_type pa = ipow<2>(pah); real_type pa2 = ipow<2>(pa); result.expnt[0] = par_[0] / (1. + par_[1] * p4 * pa) diff --git a/src/celeritas/neutron/model/NeutronElasticModel.cc b/src/celeritas/neutron/model/NeutronElasticModel.cc index 0768e514fe..c08b5b80a7 100644 --- a/src/celeritas/neutron/model/NeutronElasticModel.cc +++ b/src/celeritas/neutron/model/NeutronElasticModel.cc @@ -58,13 +58,13 @@ NeutronElasticModel::NeutronElasticModel(ActionId id, CELER_ASSERT(data.xs.elements.size() == materials.num_elements()); // Add A(Z,N)-dependent coefficients of the CHIPS elastic interaction - make_builder(&data.Q2.coeffs).reserve(materials.num_isotopes()); + make_builder(&data.diff_xs.coeffs).reserve(materials.num_isotopes()); for (auto iso_id : range(IsotopeId{materials.num_isotopes()})) { AtomicMassNumber a = materials.get(iso_id).atomic_mass_number(); - this->append_coeffs(a, &data.Q2); + this->append_coeffs(a, &data.diff_xs); } - CELER_ASSERT(data.Q2.coeffs.size() == materials.num_isotopes()); + CELER_ASSERT(data.diff_xs.coeffs.size() == materials.num_isotopes()); // Move to mirrored data, copying to device mirror_ = CollectionMirror{std::move(data)}; @@ -165,7 +165,7 @@ void NeutronElasticModel::append_xs(ImportPhysicsVector const& inp, * Wellisch, Eur. Phys. J. A9 (2001). */ void NeutronElasticModel::append_coeffs(AtomicMassNumber A, - HostDiffXsData* Q2_data) const + HostDiffXsData* data) const { ChipsDiffXsCoefficients coeffs; @@ -291,7 +291,7 @@ void NeutronElasticModel::append_coeffs(AtomicMassNumber A, coeffs.par[41] = 900 * sa / (1 + 500 / a3); } - make_builder(&Q2_data->coeffs).push_back(coeffs); + make_builder(&data->coeffs).push_back(coeffs); } //---------------------------------------------------------------------------// diff --git a/test/celeritas/neutron/NeutronElastic.test.cc b/test/celeritas/neutron/NeutronElastic.test.cc index 0895c9f185..7a1edc78d3 100644 --- a/test/celeritas/neutron/NeutronElastic.test.cc +++ b/test/celeritas/neutron/NeutronElastic.test.cc @@ -133,10 +133,10 @@ TEST_F(NeutronElasticTest, diff_xs_coeffs) { // Get A-dependent parameters of CHIPS differential cross sections used // for sampling the momentum transfer. - auto const& Q2 = model_->host_ref().Q2; + auto const& diff_xs = model_->host_ref().diff_xs; // Set the target isotope: He4 (36 parameters for light nuclei, A < 6.5) - ChipsDiffXsCoefficients he4_coeff = Q2.coeffs[IsotopeId{1}]; + ChipsDiffXsCoefficients he4_coeff = diff_xs.coeffs[IsotopeId{1}]; EXPECT_EQ(he4_coeff.par.size(), 42); EXPECT_EQ(he4_coeff.par[0], 16000); EXPECT_EQ(he4_coeff.par[10], 26.99741289550769); @@ -146,7 +146,7 @@ TEST_F(NeutronElasticTest, diff_xs_coeffs) EXPECT_EQ(he4_coeff.par[36], 0); // Set the target isotope: Cu63 (42 parameters for heavy nuclei, A > 6.5) - ChipsDiffXsCoefficients cu63_coeff = Q2.coeffs[IsotopeId{2}]; + ChipsDiffXsCoefficients cu63_coeff = diff_xs.coeffs[IsotopeId{2}]; EXPECT_EQ(cu63_coeff.par.size(), 42); EXPECT_EQ(cu63_coeff.par[0], 527.781478797624); EXPECT_EQ(cu63_coeff.par[10], 9.842761904761872); @@ -198,7 +198,7 @@ TEST_F(NeutronElasticTest, basic) TEST_F(NeutronElasticTest, extended) { - // Set the target isotope : Cu63 + // Set the target isotope : Cu63 ElementComponentId el_id{1}; IsotopeComponentId iso_id{0}; IsotopeView const isotope_he4 = this->material_track() @@ -241,9 +241,9 @@ TEST_F(NeutronElasticTest, extended) Interaction result = interact(rng_engine); // Check scattered energy and angle at each incident neutron energy - EXPECT_EQ(result.energy.value(), expected_energy[i]); - EXPECT_EQ(dot_product(result.direction, this->direction()), - expected_angle[i]); + EXPECT_SOFT_EQ(result.energy.value(), expected_energy[i]); + EXPECT_SOFT_EQ(dot_product(result.direction, this->direction()), + expected_angle[i]); energy *= factor; } }