celeritas/optical/ScintillationGenerator.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/ScintillationGenerator.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Assert.hh"
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012 #include "corecel/math/ArrayOperators.hh"
0013 #include "corecel/math/ArrayUtils.hh"
0014 #include "corecel/random/distribution/ExponentialDistribution.hh"
0015 #include "corecel/random/distribution/GenerateCanonical.hh"
0016 #include "corecel/random/distribution/NormalDistribution.hh"
0017 #include "corecel/random/distribution/RejectionSampler.hh"
0018 #include "corecel/random/distribution/Selector.hh"
0019 #include "corecel/random/distribution/UniformRealDistribution.hh"
0020
0021 #include "GeneratorDistributionData.hh"
0022 #include "ScintillationData.hh"
0023 #include "TrackInitializer.hh"
0024
0025 #include "detail/OpticalUtils.hh"
0026
0027 namespace celeritas
0028 {
0029 namespace optical
0030 {
0031 //---------------------------------------------------------------------------//
0032 /*!
0033  * Sample scintillation photons from optical property data and step data.
0034  *
0035  * The optical photons are generated evenly along the step and are emitted
0036  * uniformly over the entire solid angle with a random linear polarization.
0037  * The photon energy is calculated by the scintillation emission wavelength
0038  * \f[
0039    E = \frac{hc}{\lambda},
0040  * \f]
0041  * where \f$ h \f$ is the Planck constant and \f$ c \f$ is the speed of light,
0042  * and \f$ \lambda \f$ is sampled by the normal distribution with the mean of
0043  * scintillation emission spectrum and the standard deviation. The emitted time
0044  * is simulated according to empirical shapes of the material-dependent
0045  * scintillation time structure with one or double exponentials.
0046
0047  * \note This performs the same sampling routine as in G4Scintillation class
0048  * of the Geant4 release 11.2 with some modifications.
0049  */
0050 class ScintillationGenerator
0051 {
0052   public:
0053     //!@{
0054     //! \name Type aliases
0055     using Energy = units::MevEnergy;
0056     //!@}
0057
0058   public:
0059     // Construct from scintillation data and distribution parameters
0060     inline CELER_FUNCTION
0061     ScintillationGenerator(NativeCRef<ScintillationData> const& shared,
0062                            GeneratorDistributionData const& dist);
0063
0064     // Sample a single photon from the distribution
0065     template<class Generator>
0066     inline CELER_FUNCTION TrackInitializer operator()(Generator& rng);
0067
0068   private:
0069     //// TYPES ////
0070
0071     using UniformRealDist = UniformRealDistribution<real_type>;
0072     using ExponentialDist = ExponentialDistribution<real_type>;
0073
0074     //// DATA ////
0075
0076     GeneratorDistributionData const& dist_;
0077     NativeCRef<ScintillationData> const& shared_;
0078
0079     UniformRealDist sample_cost_;
0080     UniformRealDist sample_phi_;
0081     NormalDistribution<real_type> sample_lambda_;
0082
0083     bool is_neutral_{};
0084     units::LightSpeed delta_speed_{};
0085     Real3 delta_pos_{};
0086 };
0087
0088 //---------------------------------------------------------------------------//
0089 // INLINE DEFINITIONS
0090 //---------------------------------------------------------------------------//
0091 /*!
0092  * Construct from shared scintillation data and distribution parameters.
0093  */
0094 CELER_FUNCTION
0095 ScintillationGenerator::ScintillationGenerator(
0096     NativeCRef<ScintillationData> const& shared,
0097     GeneratorDistributionData const& dist)
0098     : dist_(dist)
0099     , shared_(shared)
0100     , sample_cost_(-1, 1)
0101     , sample_phi_(0, real_type(2 * constants::pi))
0102     , is_neutral_{dist_.charge == zero_quantity()}
0103 {
0104     if (shared_.scintillation_by_particle())
0105     {
0106         // TODO: implement sampling for particles
0107         CELER_ASSERT_UNREACHABLE();
0108     }
0109
0110     CELER_EXPECT(dist_);
0111     CELER_EXPECT(shared_);
0112
0113     auto const& pre_step = dist_.points[StepPoint::pre];
0114     auto const& post_step = dist_.points[StepPoint::post];
0115     delta_pos_ = post_step.pos - pre_step.pos;
0116     delta_speed_ = post_step.speed - pre_step.speed;
0117 }
0118
0119 //---------------------------------------------------------------------------//
0120 /*!
0121  * Sample a single scintillation photon.
0122  */
0123 template<class Generator>
0124 CELER_FUNCTION TrackInitializer ScintillationGenerator::operator()(Generator& rng)
0125 {
0126     // Sample a component
0127     ScintRecord const& component = [&] {
0128         auto const& mat = shared_.materials[dist_.material];
0129
0130         auto pdf = shared_.reals[mat.yield_pdf];
0131         auto select_idx = make_selector([&pdf](size_type i) { return pdf[i]; },
0132                                         mat.yield_pdf.size());
0133         size_type component_idx = select_idx(rng);
0134         CELER_ASSERT(component_idx < mat.components.size());
0135         return shared_.scint_records[mat.components[component_idx]];
0136     }();
0137
0138     // Sample a photon for a single scintillation component, reusing the
0139     // "spare" value that the wavelength sampler might have stored
0140     sample_lambda_
0141         = NormalDistribution{component.lambda_mean, component.lambda_sigma};
0142     ExponentialDist sample_time(real_type{1} / component.fall_time);
0143
0144     TrackInitializer photon;
0145     photon.energy = detail::wavelength_to_energy(sample_lambda_(rng));
0146
0147     // Sample direction
0148     real_type cost = sample_cost_(rng);
0149     real_type phi = sample_phi_(rng);
0150     photon.direction = from_spherical(cost, phi);
0151
0152     // Sample polarization perpendicular to the photon direction
0153     photon.polarization = [&] {
0154         Real3 temp = from_spherical(
0155             (cost > 0 ? -1 : 1) * std::sqrt(1 - ipow<2>(cost)), phi);
0156         Real3 perp = {-std::sin(phi), std::cos(phi), 0};
0157         real_type sinphi, cosphi;
0158         sincospi(UniformRealDist{}(rng), &sinphi, &cosphi);
0159         for (auto j : range(3))
0160         {
0161             temp[j] = cosphi * temp[j] + sinphi * perp[j];
0162         }
0163         return make_unit_vector(temp);
0164     }();
0165
0166     // Sample position: endpoint (collision site) if neutral, else uniform
0167     real_type u = is_neutral_ ? 1 : UniformRealDist{}(rng);
0168     photon.position = dist_.points[StepPoint::pre].pos;
0169     axpy(u, delta_pos_, &photon.position);
0170
0171     // Sample time
0172     photon.time
0173         = dist_.time
0174           + u * dist_.step_length
0175                 / (native_value_from(dist_.points[StepPoint::pre].speed)
0176                    + u * real_type(0.5) * native_value_from(delta_speed_));
0177
0178     if (component.rise_time == 0)
0179     {
0180         // Sample exponentially from fall time
0181         photon.time += sample_time(rng);
0182     }
0183     else
0184     {
0185         real_type scint_time{};
0186         real_type target;
0187         do
0188         {
0189             // Sample time exponentially by fall time, then
0190             // accept with 1 - e^{-t/rise}
0191             scint_time = sample_time(rng);
0192             target = -std::expm1(-scint_time / component.rise_time);
0193         } while (RejectionSampler(target)(rng));
0194         photon.time += scint_time;
0195     }
0196     return photon;
0197 }
0198
0199 //---------------------------------------------------------------------------//
0200 }  // namespace optical
0201 }  // namespace celeritas