added ExponentialDistribution

Framework/Random/CMakeLists.txt

@@ -8,6 +8,7 @@ set (
   CORSIKArandom_HEADERS
   RNGManager.h
   UniformRealDistribution.h
+  ExponentialDistribution.h
   )
 
 set (

Framework/Random/ExponentialDistribution.h

+/**
+ * (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

Framework/Random/UniformRealDistribution.h

-
 /**
  * (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>

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));
+  }
+}