optical/interactor/WavelengthShiftGenerator.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file celeritas/optical/interactor/WavelengthShiftGenerator.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Macros.hh"
0010 #include "corecel/Types.hh"
0011 #include "corecel/data/StackAllocator.hh"
0012 #include "corecel/random/distribution/ExponentialDistribution.hh"
0013 #include "corecel/random/distribution/UniformRealDistribution.hh"
0014 #include "geocel/random/IsotropicDistribution.hh"
0015 #include "celeritas/Types.hh"
0016 #include "celeritas/grid/NonuniformGridCalculator.hh"
0017 #include "celeritas/optical/Interaction.hh"
0018 #include "celeritas/optical/ParticleTrackView.hh"
0019 #include "celeritas/optical/TrackInitializer.hh"
0020 #include "celeritas/optical/WavelengthShiftData.hh"
0021 #include "celeritas/phys/InteractionUtils.hh"
0022
0023 namespace celeritas
0024 {
0025 namespace optical
0026 {
0027 //---------------------------------------------------------------------------//
0028 /*!
0029  * Sample optical photons from the wavelength shift process.
0030  *
0031  * A wavelength shifter absorbs an incident light and reemits secondary lights
0032  * isotropically at longer wavelengths. It usually shifts the ultraviolet
0033  * region of the radiation spectrum to the visible region, which enhances the
0034  * light collection or reduces the self-absorption of the optical production.
0035  * The number of the reemitted photons follows the Poisson distribution with
0036  * the mean number of the characteristic light production, which depends on the
0037  * optical property of wavelength shifters. The polarization of the reemitted
0038  * lights is assumed to be incoherent with respect to the polarization of the
0039  * primary optical photon.
0040  *
0041  * \note This performs the same sampling routine as in the G4OpWLS class of
0042  * the Geant4 release 11.2.
0043  */
0044 class WavelengthShiftGenerator
0045 {
0046   public:
0047     // Construct with shared and state data
0048     inline CELER_FUNCTION
0049     WavelengthShiftGenerator(NativeCRef<WavelengthShiftData> const& shared,
0050                              WlsDistributionData const& distribution);
0051
0052     // Sample a WLS photon
0053     template<class Engine>
0054     inline CELER_FUNCTION TrackInitializer operator()(Engine& rng);
0055
0056   private:
0057     //// TYPES ////
0058
0059     using Energy = units::MevEnergy;
0060
0061     //// DATA ////
0062
0063     WlsDistributionData const& distribution_;
0064     real_type time_constant_;
0065     WlsTimeProfile time_profile_;
0066     NonuniformGridCalculator calc_cdf_;
0067 };
0068
0069 //---------------------------------------------------------------------------//
0070 // INLINE DEFINITIONS
0071 //---------------------------------------------------------------------------//
0072 /*!
0073  * Construct with shared and state data.
0074  */
0075 CELER_FUNCTION
0076 WavelengthShiftGenerator::WavelengthShiftGenerator(
0077     NativeCRef<WavelengthShiftData> const& shared,
0078     WlsDistributionData const& distribution)
0079     : distribution_(distribution)
0080     , time_constant_(shared.wls_record[distribution_.material].time_constant)
0081     , time_profile_(shared.time_profile)
0082     , calc_cdf_(shared.energy_cdf[distribution_.material], shared.reals)
0083 {
0084     CELER_EXPECT(distribution_);
0085     CELER_EXPECT(distribution_.energy.value() > calc_cdf_.grid().front());
0086 }
0087
0088 //---------------------------------------------------------------------------//
0089 // INLINE DEFINITIONS
0090 //---------------------------------------------------------------------------//
0091 /*!
0092  * Sampling the wavelength shift (WLS) photons.
0093  */
0094 template<class Engine>
0095 CELER_FUNCTION TrackInitializer WavelengthShiftGenerator::operator()(Engine& rng)
0096 {
0097     // Sample wavelength shifted optical photon
0098     TrackInitializer result;
0099
0100     // Sample the emitted energy from the inverse cumulative distribution
0101     // TODO: add CDF sampler; see
0102     // https://github.com/celeritas-project/celeritas/pull/1507/files#r1844973621
0103     NonuniformGridCalculator calc_energy = calc_cdf_.make_inverse();
0104     real_type energy = calc_energy(generate_canonical(rng));
0105     if (CELER_UNLIKELY(energy > value_as<Energy>(distribution_.energy)))
0106     {
0107         // Sample a restricted energy below the incident photon energy
0108         real_type cdf_max = calc_cdf_(value_as<Energy>(distribution_.energy));
0109         UniformRealDistribution<real_type> sample_cdf(0, cdf_max);
0110         energy = calc_energy(sample_cdf(rng));
0111     }
0112     CELER_ENSURE(energy < value_as<Energy>(distribution_.energy));
0113     result.energy = Energy{energy};
0114
0115     // Use the post-step position
0116     result.position = distribution_.position;
0117
0118     // Sample the emitted photon (incoherent) direction and polarization
0119     result.direction = IsotropicDistribution()(rng);
0120     result.polarization = ExitingDirectionSampler{0, result.direction}(rng);
0121
0122     // Sample the delta time (based on the exponential relaxation)
0123     result.time
0124         = distribution_.time
0125           + (time_profile_ == WlsTimeProfile::delta
0126                  ? time_constant_
0127                  : ExponentialDistribution(real_type{1} / time_constant_)(rng));
0128
0129     CELER_ENSURE(is_soft_unit_vector(result.polarization));
0130     CELER_ENSURE(soft_zero(dot_product(result.direction, result.polarization)));
0131     return result;
0132 }
0133
0134 //---------------------------------------------------------------------------//
0135 }  // namespace optical
0136 }  // namespace celeritas