celeritas/optical/ScintillationOffload.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file celeritas/optical/ScintillationOffload.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Assert.hh"
0011 #include "corecel/Macros.hh"
0012 #include "celeritas/Quantities.hh"
0013 #include "celeritas/phys/ParticleTrackView.hh"
0014 #include "celeritas/random/distribution/NormalDistribution.hh"
0015 #include "celeritas/random/distribution/PoissonDistribution.hh"
0016 #include "celeritas/track/SimTrackView.hh"
0017
0018 #include "GeneratorDistributionData.hh"
0019 #include "OffloadData.hh"
0020 #include "ScintillationData.hh"
0021
0022 namespace celeritas
0023 {
0024 //---------------------------------------------------------------------------//
0025 /*!
0026  * Sample the number of scintillation photons to be generated.
0027  *
0028  * This populates the \c GeneratorDistributionData used by the \c
0029  * ScintillationGenerator to generate optical photons using post-step and
0030  * cached pre-step data.
0031  *
0032  * The mean number of photons is a product of the energy deposition and a
0033  * material-dependent yield fraction (photons per MeV). The actual number of
0034  * photons sampled is determined by sampling:
0035  * - for large (n > 10) mean yield, from a Gaussian distribution with a
0036  *   material-dependent spread, or
0037  * - for small yields, from a Poisson distribution.
0038  */
0039 class ScintillationOffload
0040 {
0041   public:
0042     // Construct with optical properties, scintillation, and step data
0043     inline CELER_FUNCTION
0044     ScintillationOffload(ParticleTrackView const& particle,
0045                          SimTrackView const& sim,
0046                          Real3 const& pos,
0047                          units::MevEnergy energy_deposition,
0048                          NativeCRef<optical::ScintillationData> const& shared,
0049                          OffloadPreStepData const& step_data);
0050
0051     // Populate an optical distribution data for the Scintillation Generator
0052     template<class Generator>
0053     inline CELER_FUNCTION optical::GeneratorDistributionData
0054     operator()(Generator& rng);
0055
0056   private:
0057     units::ElementaryCharge charge_;
0058     real_type step_length_;
0059     OffloadPreStepData const& pre_step_;
0060     optical::GeneratorStepData post_step_;
0061     NativeCRef<optical::ScintillationData> const& shared_;
0062     real_type mean_num_photons_{0};
0063
0064     static CELER_CONSTEXPR_FUNCTION real_type poisson_threshold()
0065     {
0066         return 10;
0067     }
0068 };
0069
0070 //---------------------------------------------------------------------------//
0071 // INLINE DEFINITIONS
0072 //---------------------------------------------------------------------------//
0073 /*!
0074  * Construct with input parameters.
0075  */
0076 CELER_FUNCTION ScintillationOffload::ScintillationOffload(
0077     ParticleTrackView const& particle,
0078     SimTrackView const& sim,
0079     Real3 const& pos,
0080     units::MevEnergy energy_deposition,
0081     NativeCRef<optical::ScintillationData> const& shared,
0082     OffloadPreStepData const& step_data)
0083     : charge_(particle.charge())
0084     , step_length_(sim.step_length())
0085     , pre_step_(step_data)
0086     , post_step_({particle.speed(), pos})
0087     , shared_(shared)
0088 {
0089     CELER_EXPECT(step_length_ > 0);
0090     CELER_EXPECT(shared_);
0091     CELER_EXPECT(pre_step_);
0092
0093     if (shared_.scintillation_by_particle())
0094     {
0095         // TODO: implement sampling for particles, assert particle data, and
0096         // cache mean number of photons
0097         CELER_ASSERT_UNREACHABLE();
0098     }
0099     else
0100     {
0101         // Scintillation will be performed on materials only
0102         CELER_ASSERT(pre_step_.material < shared_.materials.size());
0103         auto const& material = shared_.materials[pre_step_.material];
0104
0105         // TODO: Use visible energy deposition when Birks law is implemented
0106         if (material)
0107         {
0108             mean_num_photons_ = material.yield_per_energy
0109                                 * energy_deposition.value();
0110         }
0111     }
0112 }
0113
0114 //---------------------------------------------------------------------------//
0115 /*!
0116  * Return an \c GeneratorDistributionData object. If no photons are sampled, an
0117  * empty object is returned and can be verified via its own operator bool.
0118  */
0119 template<class Generator>
0120 CELER_FUNCTION optical::GeneratorDistributionData
0121 ScintillationOffload::operator()(Generator& rng)
0122 {
0123     // Material-only sampling
0124     optical::GeneratorDistributionData result;
0125     if (mean_num_photons_ > poisson_threshold())
0126     {
0127         real_type sigma = shared_.resolution_scale[pre_step_.material]
0128                           * std::sqrt(mean_num_photons_);
0129         result.num_photons = static_cast<size_type>(clamp_to_nonneg(
0130             NormalDistribution<real_type>(mean_num_photons_, sigma)(rng)
0131             + real_type{0.5}));
0132     }
0133     else if (mean_num_photons_ > 0)
0134     {
0135         result.num_photons = static_cast<size_type>(
0136             PoissonDistribution<real_type>(mean_num_photons_)(rng));
0137     }
0138
0139     if (result.num_photons > 0)
0140     {
0141         // Assign remaining data
0142         result.time = pre_step_.time;
0143         result.step_length = step_length_;
0144         result.charge = charge_;
0145         result.material = pre_step_.material;
0146         result.points[StepPoint::pre].speed = pre_step_.speed;
0147         result.points[StepPoint::pre].pos = pre_step_.pos;
0148         result.points[StepPoint::post] = post_step_;
0149     }
0150     return result;
0151 }
0152
0153 //---------------------------------------------------------------------------//
0154 }  // namespace celeritas