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