celeritas/optical/CerenkovGenerator.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/CerenkovGenerator.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011
0012 #include "corecel/Assert.hh"
0013 #include "corecel/Macros.hh"
0014 #include "corecel/Types.hh"
0015 #include "corecel/data/Collection.hh"
0016 #include "corecel/math/ArrayOperators.hh"
0017 #include "corecel/math/ArrayUtils.hh"
0018 #include "celeritas/grid/GenericCalculator.hh"
0019 #include "celeritas/random/distribution/BernoulliDistribution.hh"
0020 #include "celeritas/random/distribution/RejectionSampler.hh"
0021 #include "celeritas/random/distribution/UniformRealDistribution.hh"
0022
0023 #include "CerenkovDndxCalculator.hh"
0024 #include "GeneratorDistributionData.hh"
0025 #include "MaterialView.hh"
0026 #include "TrackInitializer.hh"
0027
0028 namespace celeritas
0029 {
0030 namespace optical
0031 {
0032 //---------------------------------------------------------------------------//
0033 /*!
0034  * Sample Cerenkov photons from the given distribution.
0035  *
0036  * Cerenkov radiation is emitted when a charged particle passes through a
0037  * dielectric medium faster than the speed of light in that medium. Photons are
0038  * emitted on the surface of a cone, with the cone angle, \f$ \theta \f$,
0039  * measured with respect to the incident particle direction. As the particle
0040  * slows down, the cone angle and the number of emitted photons decreases and
0041  * the frequency of the emitted photons increases.
0042  *
0043  * An incident charged particle with speed \f$ \beta \f$ will emit photons at
0044  * an angle \f$ \theta \f$ given by \f$ \cos\theta = 1 / (\beta n) \f$ where
0045  * \f$ n \f$ is the index of refraction of the matarial. The photon energy
0046  * \f$ \epsilon \f$ is sampled from the PDF \f[
0047    f(\epsilon) = \left[1 - \frac{1}{n^2(\epsilon)\beta^2}\right]
0048  * \f]
0049  */
0050 class CerenkovGenerator
0051 {
0052   public:
0053     // Construct from optical materials and distribution parameters
0054     inline CELER_FUNCTION
0055     CerenkovGenerator(MaterialView const& material,
0056                       NativeCRef<CerenkovData> const& shared,
0057                       GeneratorDistributionData const& dist);
0058
0059     // Sample a Cerenkov photon from the distribution
0060     template<class Generator>
0061     inline CELER_FUNCTION TrackInitializer operator()(Generator& rng);
0062
0063   private:
0064     //// TYPES ////
0065
0066     using UniformRealDist = UniformRealDistribution<real_type>;
0067
0068     //// DATA ////
0069
0070     GeneratorDistributionData const& dist_;
0071     GenericCalculator calc_refractive_index_;
0072     UniformRealDist sample_phi_;
0073     UniformRealDist sample_num_photons_;
0074     UniformRealDist sample_energy_;
0075     Real3 dir_;
0076     Real3 delta_pos_;
0077     units::LightSpeed delta_speed_;
0078     real_type delta_num_photons_;
0079     real_type dndx_pre_;
0080     real_type sin_max_sq_;
0081     real_type inv_beta_;
0082 };
0083
0084 //---------------------------------------------------------------------------//
0085 // INLINE DEFINITIONS
0086 //---------------------------------------------------------------------------//
0087 /*!
0088  * Construct from optical materials and distribution parameters.
0089  */
0090 CELER_FUNCTION
0091 CerenkovGenerator::CerenkovGenerator(MaterialView const& material,
0092                                      NativeCRef<CerenkovData> const& shared,
0093                                      GeneratorDistributionData const& dist)
0094     : dist_(dist)
0095     , calc_refractive_index_(material.make_refractive_index_calculator())
0096     , sample_phi_(0, 2 * constants::pi)
0097 {
0098     CELER_EXPECT(shared);
0099     CELER_EXPECT(dist_);
0100     CELER_EXPECT(material.material_id() == dist_.material);
0101
0102     using LS = units::LightSpeed;
0103
0104     // Calculate the mean number of photons produced per unit length at the
0105     // pre- and post-step energies
0106     auto const& pre_step = dist_.points[StepPoint::pre];
0107     auto const& post_step = dist_.points[StepPoint::post];
0108     CerenkovDndxCalculator calc_dndx(material, shared, dist_.charge);
0109     dndx_pre_ = calc_dndx(pre_step.speed);
0110     real_type dndx_post = calc_dndx(post_step.speed);
0111
0112     // Helper used to sample the displacement
0113     sample_num_photons_ = UniformRealDist(0, max(dndx_pre_, dndx_post));
0114
0115     // Helper to sample exiting photon energies
0116     auto const& energy_grid = calc_refractive_index_.grid();
0117     sample_energy_ = UniformRealDist(energy_grid.front(), energy_grid.back());
0118
0119     // Calculate 1 / beta and the max sin^2 theta
0120     inv_beta_
0121         = 2 / (value_as<LS>(pre_step.speed) + value_as<LS>(post_step.speed));
0122     CELER_ASSERT(inv_beta_ > 1);
0123     real_type cos_max = inv_beta_ / calc_refractive_index_(energy_grid.back());
0124     sin_max_sq_ = 1 - ipow<2>(cos_max);
0125
0126     // Calculate changes over the step
0127     delta_pos_ = post_step.pos - pre_step.pos;
0128     delta_num_photons_ = dndx_post - dndx_pre_;
0129     delta_speed_ = post_step.speed - pre_step.speed;
0130
0131     // Incident particle direction
0132     dir_ = make_unit_vector(delta_pos_);
0133 }
0134
0135 //---------------------------------------------------------------------------//
0136 /*!
0137  * Sample Cerenkov photons from the distribution.
0138  */
0139 template<class Generator>
0140 CELER_FUNCTION TrackInitializer CerenkovGenerator::operator()(Generator& rng)
0141 {
0142     // Sample energy and direction
0143     real_type energy;
0144     real_type cos_theta;
0145     real_type sin_theta_sq;
0146     do
0147     {
0148         // Sample an energy uniformly within the grid bounds, rejecting
0149         // if the refractive index at the sampled energy is such that the
0150         // incident particle's average speed is subluminal at that photon
0151         // wavelength.
0152         // We could improve sampling efficiency for this edge case by
0153         // increasing the minimum energy (as is done in
0154         // CerenkovDndxCalculator) to where the refractive index satisfies
0155         // this condition, but since fewer photons are emitted at lower
0156         // energies in general, relatively few rejections will take place
0157         // here.
0158         do
0159         {
0160             energy = sample_energy_(rng);
0161             cos_theta = inv_beta_ / calc_refractive_index_(energy);
0162         } while (cos_theta > 1);
0163         sin_theta_sq = 1 - ipow<2>(cos_theta);
0164     } while (RejectionSampler{sin_theta_sq, sin_max_sq_}(rng));
0165
0166     // Sample azimuthal photon direction
0167     real_type phi = sample_phi_(rng);
0168     TrackInitializer photon;
0169     photon.direction = rotate(from_spherical(cos_theta, phi), dir_);
0170     photon.energy = units::MevEnergy(energy);
0171
0172     // Photon polarization is perpendicular to the cone's surface
0173     photon.polarization
0174         = rotate(from_spherical(-std::sqrt(sin_theta_sq), phi), dir_);
0175
0176     // Sample fraction along the step
0177     UniformRealDistribution sample_step_fraction;
0178     real_type u;
0179     do
0180     {
0181         u = sample_step_fraction(rng);
0182     } while (sample_num_photons_(rng) > dndx_pre_ + u * delta_num_photons_);
0183
0184     real_type delta_time
0185         = u * dist_.step_length
0186           / (native_value_from(dist_.points[StepPoint::pre].speed)
0187              + u * real_type(0.5) * native_value_from(delta_speed_));
0188     photon.time = dist_.time + delta_time;
0189     photon.position = dist_.points[StepPoint::pre].pos;
0190     axpy(u, delta_pos_, &photon.position);
0191     return photon;
0192 }
0193
0194 //---------------------------------------------------------------------------//
0195 }  // namespace optical
0196 }  // namespace celeritas