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