Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-09-17 08:53:43

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/CherenkovDndxCalculator.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/grid/VectorUtils.hh"
0015 #include "celeritas/Quantities.hh"
0016 #include "celeritas/grid/NonuniformGridCalculator.hh"
0017 
0018 #include "CherenkovData.hh"
0019 #include "MaterialView.hh"
0020 
0021 namespace celeritas
0022 {
0023 namespace optical
0024 {
0025 //---------------------------------------------------------------------------//
0026 /*!
0027  * Calculate the mean number of Cherenkov photons produced per unit length.
0028  *
0029  * The average number of photons produced is given by
0030  * \f[
0031    \dif N = \frac{\alpha z^2}{\hbar c}\sin^2\theta \dif\epsilon \dif x =
0032    \frac{\alpha z^2}{\hbar c}\left(1 - \frac{1}{n^2\beta^2}\right) \dif\epsilon
0033    \dif x,
0034  * \f]
0035  * where \f$ n \f$ is the refractive index of the material, \f$ \epsilon \f$
0036  * is the photon energy, and \f$ \theta \f$ is the angle of the emitted photons
0037  * with respect to the incident particle direction, given by \f$ \cos\theta = 1
0038  * / (\beta n) \f$. Note that in a dispersive medium, the index of refraction
0039  * is an inreasing function of photon energy. The mean number of photons per
0040  * unit length is given by
0041  * \f[
0042    \difd{N}{x} = \frac{\alpha z^2}{\hbar c}
0043    \int_{\epsilon_\text{min}}^{\epsilon_\text{max}} \left(1 -
0044    \frac{1}{n^2\beta^2} \right) \dif\epsilon
0045    = \frac{\alpha z^2}{\hbar c}
0046    \left[\epsilon_\text{max} - \epsilon_\text{min} - \frac{1}{\beta^2}
0047    \int_{\epsilon_\text{min}}^{\epsilon_\text{max}}
0048    \frac{1}{n^2(\epsilon)}\dif\epsilon \right].
0049  * \f]
0050  */
0051 class CherenkovDndxCalculator
0052 {
0053   public:
0054     // Construct from optical materials and Cherenkov angle integrals
0055     inline CELER_FUNCTION
0056     CherenkovDndxCalculator(MaterialView const& material,
0057                            NativeCRef<CherenkovData> const& shared,
0058                            units::ElementaryCharge charge);
0059 
0060     // Calculate the mean number of Cherenkov photons produced per unit length
0061     inline CELER_FUNCTION real_type operator()(units::LightSpeed beta);
0062 
0063   private:
0064     // Calculate refractive index [MeV -> unitless]
0065     NonuniformGridCalculator calc_refractive_index_;
0066 
0067     // Calculate the Cherenkov angle integral [MeV -> unitless]
0068     NonuniformGridCalculator calc_integral_;
0069 
0070     // Square of particle charge
0071     real_type zsq_;
0072 };
0073 
0074 //---------------------------------------------------------------------------//
0075 // INLINE DEFINITIONS
0076 //---------------------------------------------------------------------------//
0077 /*!
0078  * Construct from optical materials and Cherenkov angle integrals.
0079  */
0080 CELER_FUNCTION
0081 CherenkovDndxCalculator::CherenkovDndxCalculator(
0082     MaterialView const& material,
0083     NativeCRef<CherenkovData> const& shared,
0084     units::ElementaryCharge charge)
0085     : calc_refractive_index_(material.make_refractive_index_calculator())
0086     , calc_integral_{shared.angle_integral[material.material_id()],
0087                      shared.reals}
0088     , zsq_(ipow<2>(charge.value()))
0089 {
0090     CELER_EXPECT(charge != zero_quantity());
0091 }
0092 
0093 //---------------------------------------------------------------------------//
0094 /*!
0095  * Calculate the mean number of Cherenkov photons produced per unit length.
0096  *
0097  * \todo define a "generic grid extents" class for finding the lower/upper x/y
0098  * coordinates on grid data. In the future we could cache these if the memory
0099  * lookups result in too much indirection.
0100  */
0101 CELER_FUNCTION real_type
0102 CherenkovDndxCalculator::operator()(units::LightSpeed beta)
0103 {
0104     CELER_EXPECT(beta.value() > 0 && beta.value() <= 1);
0105 
0106     real_type inv_beta = 1 / beta.value();
0107     real_type energy_max = calc_refractive_index_.grid().back();
0108     if (inv_beta > calc_refractive_index_(energy_max))
0109     {
0110         // Incident particle energy is below the threshold for Cherenkov
0111         // emission (i.e., beta < 1 / n_max)
0112         return 0;
0113     }
0114 
0115     // Calculate \f$ \int_{\epsilon_\text{min}}^{\epsilon_\text{max}}
0116     // \dif\epsilon \left(1 - \frac{1}{n^2\beta^2}\right) \f$
0117     real_type energy;
0118     if (inv_beta < calc_refractive_index_[0])
0119     {
0120         energy = energy_max - calc_refractive_index_.grid().front()
0121                  - calc_integral_(energy_max) * ipow<2>(inv_beta);
0122     }
0123     else
0124     {
0125         // Find the energy where the refractive index is equal to 1 / beta.
0126         // Both energy and refractive index are monotonically increasing, so
0127         // the grid and values can be swapped and the energy can be calculated
0128         // from a given index of refraction
0129         real_type energy_min = calc_refractive_index_.make_inverse()(inv_beta);
0130         energy = energy_max - energy_min
0131                  - (calc_integral_(energy_max) - calc_integral_(energy_min))
0132                        * ipow<2>(inv_beta);
0133     }
0134 
0135     // Calculate number of photons. This may be negative if the incident
0136     // particle energy is very close to (just above) the Cherenkov production
0137     // threshold
0138     return clamp_to_nonneg(zsq_
0139                            * (constants::alpha_fine_structure
0140                               / (constants::hbar_planck * constants::c_light))
0141                            * native_value_from(units::MevEnergy(energy)));
0142 }
0143 
0144 //---------------------------------------------------------------------------//
0145 }  // namespace optical
0146 }  // namespace celeritas