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