optical/gen/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/gen/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 "GeneratorData.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<ScintillationData> const& shared,
0048                          OffloadPreStepData const& step_data);
0049
0050     // Gather the input data needed to sample scintillation photons
0051     template<class Generator>
0052     inline CELER_FUNCTION GeneratorDistributionData operator()(Generator& rng);
0053
0054   private:
0055     units::ElementaryCharge charge_;
0056     real_type step_length_;
0057     OffloadPreStepData const& pre_step_;
0058     GeneratorStepData post_step_;
0059     NativeCRef<ScintillationData> const& shared_;
0060     real_type mean_num_photons_{0};
0061
0062     static CELER_CONSTEXPR_FUNCTION real_type poisson_threshold()
0063     {
0064         return 10;
0065     }
0066 };
0067
0068 //---------------------------------------------------------------------------//
0069 // INLINE DEFINITIONS
0070 //---------------------------------------------------------------------------//
0071 /*!
0072  * Construct with input parameters.
0073  */
0074 CELER_FUNCTION ScintillationOffload::ScintillationOffload(
0075     ParticleTrackView const& particle,
0076     SimTrackView const& sim,
0077     Real3 const& pos,
0078     units::MevEnergy energy_deposition,
0079     NativeCRef<ScintillationData> const& shared,
0080     OffloadPreStepData const& step_data)
0081     : charge_(particle.charge())
0082     , step_length_(sim.step_length())
0083     , pre_step_(step_data)
0084     , post_step_({particle.speed(), pos})
0085     , shared_(shared)
0086 {
0087     CELER_EXPECT(step_length_ > 0);
0088     CELER_EXPECT(shared_);
0089     CELER_EXPECT(pre_step_);
0090
0091     if (shared_.scintillation_by_particle())
0092     {
0093         //! \todo Implement sampling for particles
0094         CELER_ASSERT_UNREACHABLE();
0095     }
0096     else
0097     {
0098         // Scintillation will be performed on materials only
0099         CELER_ASSERT(pre_step_.material < shared_.materials.size());
0100         auto const& material = shared_.materials[pre_step_.material];
0101
0102         //! \todo Use visible energy deposition when Birks law is implemented
0103         if (material)
0104         {
0105             mean_num_photons_ = material.yield_per_energy
0106                                 * energy_deposition.value();
0107         }
0108     }
0109 }
0110
0111 //---------------------------------------------------------------------------//
0112 /*!
0113  * Collect the distribution data needed to sample scintillation photons.
0114  *
0115  * If no photons are sampled an empty object is returned.
0116  */
0117 template<class Generator>
0118 CELER_FUNCTION GeneratorDistributionData
0119 ScintillationOffload::operator()(Generator& rng)
0120 {
0121     // Material-only sampling
0122     GeneratorDistributionData result;
0123     if (mean_num_photons_ > poisson_threshold())
0124     {
0125         real_type sigma = shared_.resolution_scale[pre_step_.material]
0126                           * std::sqrt(mean_num_photons_);
0127         result.num_photons = static_cast<size_type>(clamp_to_nonneg(
0128             NormalDistribution<real_type>(mean_num_photons_, sigma)(rng)
0129             + real_type{0.5}));
0130     }
0131     else if (mean_num_photons_ > 0)
0132     {
0133         result.num_photons = static_cast<size_type>(
0134             PoissonDistribution<real_type>(mean_num_photons_)(rng));
0135     }
0136
0137     if (result.num_photons > 0)
0138     {
0139         // Assign remaining data
0140         result.time = pre_step_.time;
0141         result.step_length = step_length_;
0142         result.charge = charge_;
0143         result.material = pre_step_.material;
0144         result.points[StepPoint::pre].speed = pre_step_.speed;
0145         result.points[StepPoint::pre].pos = pre_step_.pos;
0146         result.points[StepPoint::post] = post_step_;
0147     }
0148     return result;
0149 }
0150
0151 //---------------------------------------------------------------------------//
0152 }  // namespace celeritas