em/distribution/EnergyLossGammaDistribution.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2022-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/em/distribution/EnergyLossGammaDistribution.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/math/Algorithms.hh"
0016 #include "celeritas/random/distribution/GammaDistribution.hh"
0017
0018 #include "EnergyLossHelper.hh"
0019
0020 namespace celeritas
0021 {
0022 //---------------------------------------------------------------------------//
0023 /*!
0024  * Sample energy loss from a gamma distribution.
0025  *
0026  * This model is valid for heavy particles with small mean losses (less than 2
0027  * Bohr standard deviations). It is a special case of
0028  * \c EnergyLossGaussianDistribution (see that class for more documentation).
0029  *
0030  * Note that while this appears in G4UniversalFluctuation, the Geant4
0031  * documentation does not explain why the loss is sampled from a gamma
0032  * distribution in this case.
0033  */
0034 class EnergyLossGammaDistribution
0035 {
0036   public:
0037     //!@{
0038     //! \name Type aliases
0039     using Energy = units::MevEnergy;
0040     using EnergySq = Quantity<UnitProduct<units::Mev, units::Mev>>;
0041     //!@}
0042
0043   public:
0044     // Construct from distribution parameters
0045     inline CELER_FUNCTION
0046     EnergyLossGammaDistribution(Energy mean_loss, EnergySq bohr_var);
0047
0048     // Construct from helper-calculated data
0049     explicit inline CELER_FUNCTION
0050     EnergyLossGammaDistribution(EnergyLossHelper const& helper);
0051
0052     //! Sample energy loss according to the distribution
0053     template<class Generator>
0054     CELER_FUNCTION Energy operator()(Generator& rng)
0055     {
0056         return Energy{sample_gamma_(rng)};
0057     }
0058
0059   private:
0060     using GammaDist = GammaDistribution<real_type>;
0061
0062     GammaDist sample_gamma_;
0063
0064     static inline CELER_FUNCTION GammaDist build_gamma(real_type mean,
0065                                                        real_type stddev);
0066 };
0067
0068 //---------------------------------------------------------------------------//
0069 // INLINE DEFINITIONS
0070 //---------------------------------------------------------------------------//
0071 /*!
0072  * Construct from distribution parameters.
0073  *
0074  * The model is only valid (rather, used in Geant4) for mean < 2 * stddev, but
0075  * it is allowable to construct the sampler explicitly but outside this range,
0076  * for analysis purposes.
0077  */
0078 CELER_FUNCTION
0079 EnergyLossGammaDistribution::EnergyLossGammaDistribution(Energy mean_loss,
0080                                                          EnergySq bohr_var)
0081     : sample_gamma_(EnergyLossGammaDistribution::build_gamma(mean_loss.value(),
0082                                                              bohr_var.value()))
0083 {
0084     CELER_EXPECT(mean_loss > zero_quantity());
0085     CELER_EXPECT(bohr_var > zero_quantity());
0086 }
0087
0088 //---------------------------------------------------------------------------//
0089 /*!
0090  * Construct from helper-calculated data.
0091  */
0092 CELER_FUNCTION EnergyLossGammaDistribution::EnergyLossGammaDistribution(
0093     EnergyLossHelper const& helper)
0094     : EnergyLossGammaDistribution(helper.mean_loss(), helper.bohr_variance())
0095 {
0096     CELER_ASSERT(helper.model() == EnergyLossFluctuationModel::gamma);
0097 }
0098
0099 //---------------------------------------------------------------------------//
0100 /*!
0101  * Helper function to construct gamma distribution.
0102  */
0103 CELER_FUNCTION auto
0104 EnergyLossGammaDistribution::build_gamma(real_type mean,
0105                                          real_type var) -> GammaDist
0106 {
0107     real_type k = ipow<2>(mean) / var;
0108     return GammaDist{k, mean / k};
0109 }
0110
0111 //---------------------------------------------------------------------------//
0112 }  // namespace celeritas