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