diff --git a/Framework/Random/CMakeLists.txt b/Framework/Random/CMakeLists.txt index 3607b585a62dcc418c986852428fc8fa8b518325..84479e1230a92ef5a1773ff75ebdf72347c3f7db 100644 --- a/Framework/Random/CMakeLists.txt +++ b/Framework/Random/CMakeLists.txt @@ -8,6 +8,7 @@ set ( CORSIKArandom_HEADERS RNGManager.h UniformRealDistribution.h + ExponentialDistribution.h ) set ( diff --git a/Framework/Random/ExponentialDistribution.h b/Framework/Random/ExponentialDistribution.h new file mode 100644 index 0000000000000000000000000000000000000000..248d785a748651bdd3cd462e2541bbd07403d8d0 --- /dev/null +++ b/Framework/Random/ExponentialDistribution.h @@ -0,0 +1,38 @@ +/** + * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + * + * See file AUTHORS for a list of contributors. + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#ifndef _include_ExponentialDistribution_h +#define _include_ExponentialDistribution_h + +#include <corsika/units/PhysicalUnits.h> +#include <random> + +namespace corsika::random { + + template <class TQuantity> + class ExponentialDistribution { + using RealType = typename TQuantity::value_type; + std::exponential_distribution<RealType> dist{1.}; + + TQuantity const fBeta; + + public: + ExponentialDistribution(TQuantity beta) + : fBeta(beta) {} + + template <class Generator> + TQuantity operator()(Generator& g) { + return fBeta * dist(g); + } + }; + +} // namespace corsika::random + +#endif diff --git a/Framework/Random/UniformRealDistribution.h b/Framework/Random/UniformRealDistribution.h index 738729282ff41a9a3c7a5415c7a99eae93e23e38..fd501ab36b2d0e411a1a2fc256734f5ec0221932 100644 --- a/Framework/Random/UniformRealDistribution.h +++ b/Framework/Random/UniformRealDistribution.h @@ -1,4 +1,3 @@ - /** * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu * @@ -9,8 +8,8 @@ * the license. */ -#ifndef _include_random_distributions_h -#define _include_random_distributions_h +#ifndef _include_UniformRealDistribution_h +#define _include_UniformRealDistribution_h #include <corsika/units/PhysicalUnits.h> #include <random> diff --git a/Framework/Random/testRandom.cc b/Framework/Random/testRandom.cc index 56d8157a74530d404f5d8ba1c8cf8aefe9f75f99..97923eb09e19a7e6bc1cd8a4c8b6971907af9528 100644 --- a/Framework/Random/testRandom.cc +++ b/Framework/Random/testRandom.cc @@ -14,10 +14,12 @@ #include <corsika/random/RNGManager.h> #include <corsika/random/UniformRealDistribution.h> +#include <corsika/random/ExponentialDistribution.h> #include <corsika/units/PhysicalUnits.h> #include <iostream> #include <limits> #include <random> +#include <type_traits> using namespace corsika::random; @@ -83,3 +85,25 @@ TEST_CASE("UniformRealDistribution") { CHECK(max / 18_cm == Approx(1.)); } } + +TEST_CASE("ExponentialDistribution") { + using namespace corsika::units::si; + std::mt19937 rng; + + auto const beta = 15_m; + + corsika::random::ExponentialDistribution dist(beta); + + SECTION("mean") { + std::remove_const<decltype(beta)>::type mean = beta * 0; + + int constexpr N = 1'000'000; + + for (int i{0}; i < N; ++i) { + decltype(beta) x = dist(rng); + mean += x / N; + } + + CHECK(mean / beta == Approx(1).margin(1e-2)); + } +}