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