optical/interactor/WavelengthShiftInteractor.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/WavelengthShiftInteractor.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/PoissonDistribution.hh"
0014 #include "corecel/random/distribution/UniformRealDistribution.hh"
0015 #include "geocel/random/IsotropicDistribution.hh"
0016 #include "celeritas/Types.hh"
0017 #include "celeritas/grid/NonuniformGridCalculator.hh"
0018 #include "celeritas/optical/Interaction.hh"
0019 #include "celeritas/optical/ParticleTrackView.hh"
0020 #include "celeritas/optical/TrackInitializer.hh"
0021 #include "celeritas/optical/WavelengthShiftData.hh"
0022 #include "celeritas/phys/InteractionUtils.hh"
0023
0024 namespace celeritas
0025 {
0026 namespace optical
0027 {
0028 //---------------------------------------------------------------------------//
0029 /*!
0030  * Apply the wavelength shift (WLS) to optical photons.
0031  *
0032  * A wavelength shifter absorbs an incident light and reemits secondary lights
0033  * isotropically at longer wavelengths. It usually shifts the ultraviolet
0034  * region of the radiation spectrum to the visible region, which enhances the
0035  * light collection or reduces the self-absorption of the optical production.
0036  * The number of the reemitted photons follows the Poisson distribution with
0037  * the mean number of the characteristic light production, which depends on the
0038  * optical property of wavelength shifters. The polarization of the reemitted
0039  * lights is assumed to be incoherent with respect to the polarization of the
0040  * primary optical photon.
0041  *
0042  * \note This performs the same sampling routine as in the G4OpWLS class of
0043  * the Geant4 release 11.2.
0044  */
0045 class WavelengthShiftInteractor
0046 {
0047   public:
0048     //!@{
0049     //! \name Type aliases
0050     using Energy = units::MevEnergy;
0051     using ParamsRef = NativeCRef<WavelengthShiftData>;
0052     using SecondaryAllocator = StackAllocator<TrackInitializer>;
0053     //!@}
0054
0055   public:
0056     // Construct with shared and state data
0057     inline CELER_FUNCTION
0058     WavelengthShiftInteractor(ParamsRef const& shared,
0059                               ParticleTrackView const& particle,
0060                               OptMatId const& mat_id,
0061                               SecondaryAllocator& allocate);
0062
0063     // Sample an interaction with the given RNG
0064     template<class Engine>
0065     inline CELER_FUNCTION Interaction operator()(Engine& rng);
0066
0067   private:
0068     //// DATA ////
0069
0070     // Incident photon energy
0071     Energy const inc_energy_;
0072     // Sample distributions
0073     PoissonDistribution<real_type> sample_num_photons_;
0074     ExponentialDistribution<real_type> sample_time_;
0075     // Grid calculators
0076     NonuniformGridCalculator calc_cdf_;
0077     // Allocate space for secondary particles
0078     SecondaryAllocator& allocate_;
0079 };
0080
0081 //---------------------------------------------------------------------------//
0082 // INLINE DEFINITIONS
0083 //---------------------------------------------------------------------------//
0084 /*!
0085  * Construct with shared and state data.
0086  */
0087 CELER_FUNCTION
0088 WavelengthShiftInteractor::WavelengthShiftInteractor(
0089     ParamsRef const& shared,
0090     ParticleTrackView const& particle,
0091     OptMatId const& mat_id,
0092     SecondaryAllocator& allocate)
0093     : inc_energy_(particle.energy())
0094     , sample_num_photons_(shared.wls_record[mat_id].mean_num_photons)
0095     , sample_time_(real_type{1} / shared.wls_record[mat_id].time_constant)
0096     , calc_cdf_(shared.energy_cdf[mat_id], shared.reals)
0097     , allocate_(allocate)
0098 {
0099     CELER_EXPECT(inc_energy_.value() > calc_cdf_.grid().front());
0100 }
0101
0102 //---------------------------------------------------------------------------//
0103 // INLINE DEFINITIONS
0104 //---------------------------------------------------------------------------//
0105 /*!
0106  * Sampling the wavelength shift (WLS) photons.
0107  */
0108 template<class Engine>
0109 CELER_FUNCTION Interaction WavelengthShiftInteractor::operator()(Engine& rng)
0110 {
0111     /*!
0112      * Sample the number of photons generated from WLS.
0113      *
0114      * \todo if this is nonzero and allocation fails, we lose reproducibility.
0115      */
0116     unsigned int num_photons = sample_num_photons_(rng);
0117
0118     if (num_photons == 0)
0119     {
0120         // Return absorbed photon without reemitted optical photons
0121         return Interaction::from_absorption();
0122     }
0123
0124     // Allocate space for reemitted optical photons - Note: the reproducibility
0125     // is not guaranteed in the case of an out-of-memory error
0126     TrackInitializer* secondaries = allocate_(num_photons);
0127
0128     if (secondaries == nullptr)
0129     {
0130         // Failed to allocate space
0131         return Interaction::from_failure();
0132     }
0133
0134     // Sample wavelength shifted optical photons
0135     Interaction result = Interaction::from_absorption();
0136     result.secondaries = {secondaries, num_photons};
0137
0138     IsotropicDistribution sample_direction{};
0139     NonuniformGridCalculator calc_energy = calc_cdf_.make_inverse();
0140     for (size_type i : range(num_photons))
0141     {
0142         // Sample the emitted energy from the inverse cumulative distribution
0143         // TODO: add CDF sampler; see
0144         // https://github.com/celeritas-project/celeritas/pull/1507/files#r1844973621
0145         real_type energy = calc_energy(generate_canonical(rng));
0146         if (CELER_UNLIKELY(energy > inc_energy_.value()))
0147         {
0148             // Sample a restricted energy below the incident photon energy
0149             real_type cdf_max = calc_cdf_(inc_energy_.value());
0150             UniformRealDistribution<real_type> sample_cdf(0, cdf_max);
0151             energy = calc_energy(sample_cdf(rng));
0152         }
0153         CELER_ENSURE(energy < inc_energy_.value());
0154         secondaries[i].energy = Energy{energy};
0155
0156         // Sample the emitted photon (incoherent) direction and polarization
0157         secondaries[i].direction = sample_direction(rng);
0158         secondaries[i].polarization
0159             = ExitingDirectionSampler{0, secondaries[i].direction}(rng);
0160
0161         CELER_ENSURE(is_soft_unit_vector(secondaries[i].direction));
0162         CELER_ENSURE(is_soft_unit_vector(secondaries[i].polarization));
0163         CELER_ENSURE(soft_zero(dot_product(secondaries[i].direction,
0164                                            secondaries[i].polarization)));
0165
0166         // Sample the delta time (based on the exponential relaxation)
0167         secondaries[i].time = sample_time_(rng);
0168     }
0169
0170     CELER_ENSURE(result.action == Interaction::Action::absorbed);
0171
0172     return result;
0173 }
0174
0175 //---------------------------------------------------------------------------//
0176 }  // namespace optical
0177 }  // namespace celeritas