Skip to content

Commit

Permalink
Add Wentzel OK&VI transport cross section calculator (celeritas-proje…
Browse files Browse the repository at this point in the history 
…ct#1116)

* Use MevEnergy for Wentzel incident particle cutoff energy
* Add Wentzel OK and VI transport xs calculator
* Add Wentzel transport xs calculator test
* Rename WentzelRatioCalculator --> WentzelHelper
* Address review feedback
  • Loading branch information
@amandalund
amandalund committed Feb 20, 2024
1 parent a8e1ac4 commit 4d47e34
Show file tree
Hide file tree
Showing 5 changed files with 331 additions and 99 deletions.
18 changes: 9 additions & 9 deletions src/celeritas/em/distribution/WentzelDistribution.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
#include "celeritas/em/data/WentzelData.hh"
#include "celeritas/em/interactor/detail/PhysicsConstants.hh"
#include "celeritas/em/xs/MottRatioCalculator.hh"
#include "celeritas/em/xs/WentzelRatioCalculator.hh"
#include "celeritas/em/xs/WentzelHelper.hh"
#include "celeritas/mat/IsotopeView.hh"
#include "celeritas/phys/ParticleTrackView.hh"
#include "celeritas/random/distribution/BernoulliDistribution.hh"
Expand Down Expand Up @@ -56,7 +56,7 @@ class WentzelDistribution
WentzelDistribution(ParticleTrackView const& particle,
IsotopeView const& target,
WentzelElementData const& element_data,
real_type cutoff_energy,
Energy cutoff,
WentzelRef const& data);

// Sample the polar scattering angle
Expand All @@ -78,8 +78,8 @@ class WentzelDistribution
// Mott coefficients for the target element
WentzelElementData const& element_data_;

// Ratio of elecron to total cross sections for the Wentzel model
WentzelRatioCalculator const calc_elec_ratio_;
// Helper for calculating xs ratio and other quantities
WentzelHelper const helper_;

// Calculates the form factor from the scattered polar angle
inline CELER_FUNCTION real_type calculate_form_factor(real_type cos_t) const;
Expand Down Expand Up @@ -115,13 +115,13 @@ CELER_FUNCTION
WentzelDistribution::WentzelDistribution(ParticleTrackView const& particle,
IsotopeView const& target,
WentzelElementData const& element_data,
real_type cutoff_energy,
Energy cutoff,
WentzelRef const& data)
: data_(data)
, particle_(particle)
, target_(target)
, element_data_(element_data)
, calc_elec_ratio_(particle, target.atomic_number(), data, cutoff_energy)
, helper_(particle, target.atomic_number(), data, cutoff)
{
}

Expand All @@ -133,10 +133,10 @@ template<class Engine>
CELER_FUNCTION real_type WentzelDistribution::operator()(Engine& rng) const
{
real_type cos_theta = 1;
if (BernoulliDistribution(calc_elec_ratio_())(rng))
if (BernoulliDistribution(helper_.calc_xs_ratio())(rng))
{
// Scattered off of electrons
cos_theta = this->sample_cos_t(calc_elec_ratio_.cos_t_max_elec(), rng);
cos_theta = this->sample_cos_t(helper_.costheta_max_electron(), rng);
}
else
{
Expand Down Expand Up @@ -259,7 +259,7 @@ CELER_FUNCTION real_type WentzelDistribution::sample_cos_t(real_type cos_t_max,
// 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 sc = calc_elec_ratio_.screening_coefficient();
real_type const sc = helper_.screening_coefficient();

return clamp(1 + 2 * sc * mu * (1 - xi) / (sc - mu * xi),
real_type{-1},
Expand Down
2 changes: 1 addition & 1 deletion src/celeritas/em/interactor/WentzelInteractor.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,7 @@ WentzelInteractor::WentzelInteractor(WentzelRef const& shared,
target,
shared.elem_data[el_id],
// TODO: Use proton when supported
value_as<Energy>(cutoffs.energy(shared.ids.electron)),
cutoffs.energy(shared.ids.electron),
shared)
{
CELER_EXPECT(particle_.particle_id() == shared.ids.electron
Expand Down
152 changes: 78 additions & 74 deletions ...celeritas/em/xs/WentzelRatioCalculator.hh → src/celeritas/em/xs/WentzelHelper.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: (Apache-2.0 OR MIT)
//---------------------------------------------------------------------------//
//! \file celeritas/em/xs/WentzelRatioCalculator.hh
//! \file celeritas/em/xs/WentzelHelper.hh
//---------------------------------------------------------------------------//
#pragma once

Expand All @@ -17,10 +17,16 @@ namespace celeritas
{
//---------------------------------------------------------------------------//
/*!
* Calculate the ratio of the electron to total Wentzel cross sections for
* elastic Coulomb scattering.
* Helper class for the Wentzel OK and VI Coulomb scattering model.
*
* This calculates the Moliere screening coefficient, the maximum scattering
* angle off of electrons, and the ratio of the electron to total Wentzel cross
* sections.
*
* References:
* [PRM] Geant4 Physics Reference Manual (Release 11.1) section 8.5.
*/
class WentzelRatioCalculator
class WentzelHelper
{
public:
//!@{
Expand All @@ -31,23 +37,33 @@ class WentzelRatioCalculator
//!@}

public:
//! Construct the calculator from the given values
inline CELER_FUNCTION
WentzelRatioCalculator(ParticleTrackView const& particle,
AtomicNumber target_z,
WentzelRef const& data,
real_type cutoff_energy);
// Construct from particle and material properties
inline CELER_FUNCTION WentzelHelper(ParticleTrackView const& particle,
AtomicNumber target_z,
WentzelRef const& data,
Energy cutoff);

// The ratio of electron to total cross section for Coulomb scattering
inline CELER_FUNCTION real_type operator()() const;
//! Get the target atomic number
CELER_FUNCTION AtomicNumber atomic_number() const { return target_z_; }

//! Get the Moliere screening coefficient
CELER_FUNCTION real_type screening_coefficient() const
{
return screening_coefficient_;
}

// Moilere screening coefficient
inline CELER_FUNCTION real_type screening_coefficient() const;
//! Get the maximum scattering angle off of electrons
CELER_FUNCTION real_type costheta_max_electron() const
{
return cos_t_max_elec_;
}

// (Cosine of) the maximum scattering angle off of electrons
inline CELER_FUNCTION real_type cos_t_max_elec() const;
// The ratio of electron to total cross section for Coulomb scattering
inline CELER_FUNCTION real_type calc_xs_ratio() const;

private:
//// DATA ////

// Target atomic number
AtomicNumber const target_z_;

Expand All @@ -57,49 +73,45 @@ class WentzelRatioCalculator
// Cosine of the maximum scattering angle off of electrons
real_type cos_t_max_elec_;

// Shared WentzelModel data
WentzelRef const& data_;
//// HELPER FUNCTIONS ////

//! Calculate the Moilere screening coefficient
// Calculate the Moliere screening coefficient
inline CELER_FUNCTION real_type
calc_screening_coefficient(ParticleTrackView const& particle) const;

//! Calculate the screening coefficient R^2 for electrons
// Calculate the screening coefficient R^2 for electrons
CELER_CONSTEXPR_FUNCTION real_type screen_r_sq_elec() const;

//! Calculate the (cosine of) the maximum scattering angle off of electrons
inline CELER_FUNCTION real_type calc_max_electron_cos_t(
ParticleTrackView const& particle, real_type cutoff_energy) const;
// Calculate the (cosine of) the maximum scattering angle off of electrons
inline CELER_FUNCTION real_type calc_costheta_max_electron(
ParticleTrackView const&, WentzelRef const&, Energy) const;
};

//---------------------------------------------------------------------------//
// INLINE DEFINITIONS
//---------------------------------------------------------------------------//
/*!
* Construct with state data.
* Construct from particle and material properties.
*/
CELER_FUNCTION
WentzelRatioCalculator::WentzelRatioCalculator(ParticleTrackView const& particle,
AtomicNumber target_z,
WentzelRef const& data,
real_type cutoff_energy)
: target_z_(target_z), data_(data)
WentzelHelper::WentzelHelper(ParticleTrackView const& particle,
AtomicNumber target_z,
WentzelRef const& data,
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))
{
screening_coefficient_ = calc_screening_coefficient(particle)
* data_.screening_factor;
cos_t_max_elec_ = calc_max_electron_cos_t(particle, cutoff_energy);

CELER_EXPECT(target_z_.get() > 0);
CELER_EXPECT(screening_coefficient_ > 0);
CELER_EXPECT(cos_t_max_elec_ >= -1 && cos_t_max_elec_ <= 1);
}

//---------------------------------------------------------------------------//
/*!
* Ratio of electron cross section to the total (nuclear + electron)
* cross section.
* Ratio of electron cross section to total (nuclear + electron) cross section.
*/
CELER_FUNCTION real_type WentzelRatioCalculator::operator()() const
CELER_FUNCTION real_type WentzelHelper::calc_xs_ratio() const
{
// Calculating only reduced cross sections by elimination mutual factors
// in the ratio.
Expand All @@ -112,18 +124,9 @@ CELER_FUNCTION real_type WentzelRatioCalculator::operator()() const

//---------------------------------------------------------------------------//
/*!
* Retrieve the cached Moilere screening coefficient.
* Calculate the Moliere screening coefficient as in [PRM] eqn 8.51.
*/
CELER_FUNCTION real_type WentzelRatioCalculator::screening_coefficient() const
{
return screening_coefficient_;
}

//---------------------------------------------------------------------------//
/*!
* Calculate the Moilere screening coefficient as in [PRM] eqn 8.51.
*/
CELER_FUNCTION real_type WentzelRatioCalculator::calc_screening_coefficient(
CELER_FUNCTION real_type WentzelHelper::calc_screening_coefficient(
ParticleTrackView const& particle) const
{
// TODO: Reference for just proton correction?
Expand All @@ -144,19 +147,28 @@ CELER_FUNCTION real_type WentzelRatioCalculator::calc_screening_coefficient(
* factor / particle.beta_sq());
}

return correction * screen_r_sq_elec() * sq_cbrt_z
return correction * this->screen_r_sq_elec() * sq_cbrt_z
/ value_as<MomentumSq>(particle.momentum_sq());
}

//---------------------------------------------------------------------------//
/*!
* Calculate the screening R^2 coefficient for incident electrons. This is
* the constant prefactor of [PRM] eqn 8.51.
* Calculate the constant factor of the screening coefficient.
*
* This is the constant prefactor \f$ R^2 / Z^{2/3} \f$ of the screening
* coefficient for incident electrons (equation 8.51 in [PRM]). The screening
* radius \f$ R \f$ is given by:
* \f[
R = \frac{\hbar Z^{1/3}}{2C_{TF} a_0},
* \f]
* where the Thomas-Fermi constant \f$ C_{TF} \f$ is defined as
* \f[
C_{TF} = \frac{1}{2} \left(\frac{3\pi}{4}\right)^{2/3}.
* \f]
*/
CELER_CONSTEXPR_FUNCTION real_type WentzelRatioCalculator::screen_r_sq_elec() const
CELER_CONSTEXPR_FUNCTION real_type WentzelHelper::screen_r_sq_elec() const
{
//! Thomas-Fermi constant C_TF
//! \f$ \frac{1}{2}\left(\frac{3\pi}{4}\right)^{2/3} \f$
//! Thomas-Fermi constant \f$ C_{TF} \f$
constexpr real_type ctf = 0.8853413770001135;

return native_value_to<MomentumSq>(
Expand All @@ -166,29 +178,24 @@ CELER_CONSTEXPR_FUNCTION real_type WentzelRatioCalculator::screen_r_sq_elec() co

//---------------------------------------------------------------------------//
/*!
* (Cosine of) the maximum polar angle that incident particles scatter off
* of the target's electrons.
* Calculate the maximum scattering angle off the target's electrons.
*
* 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 WentzelRatioCalculator::cos_t_max_elec() const
{
return cos_t_max_elec_;
}

//---------------------------------------------------------------------------//
/*!
* Calculate the (cosine of the) maximum polar angle that incident particle
* can scatter off of the target's electrons.
*/
CELER_FUNCTION real_type WentzelRatioCalculator::calc_max_electron_cos_t(
ParticleTrackView const& particle, real_type cutoff_energy) const
CELER_FUNCTION real_type
WentzelHelper::calc_costheta_max_electron(ParticleTrackView const& particle,
WentzelRef const& data,
Energy cutoff) const
{
real_type inc_energy = value_as<Energy>(particle.energy());
real_type mass = value_as<Mass>(particle.mass());

real_type max_energy = (particle.particle_id() == data_.ids.electron)
real_type max_energy = particle.particle_id() == data.ids.electron
? real_type{0.5} * inc_energy
: inc_energy;
real_type final_energy = inc_energy - min(cutoff_energy, max_energy);
real_type final_energy = inc_energy
- min(value_as<Energy>(cutoff), max_energy);

if (final_energy > 0)
{
Expand All @@ -198,10 +205,7 @@ CELER_FUNCTION real_type WentzelRatioCalculator::calc_max_electron_cos_t(

return clamp(cos_t_max, real_type{0}, real_type{1});
}
else
{
return 0;
}
return 0;
}

//---------------------------------------------------------------------------//
Expand Down

0 comments on commit 4d47e34

Please sign in to comment.