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