VecGeom/base/RNG.h

0001 /// \file RNG.h
0002 /// \author Johannes de Fine Licht (johannes.definelicht@cern.ch)
0003
0004 #ifndef VECGEOM_BASE_RNG_H_
0005 #define VECGEOM_BASE_RNG_H_
0006
0007 #include "VecGeom/base/Global.h"
0008
0009 #include <random>
0010
0011 namespace vecgeom {
0012
0013 /**
0014  * @brief Singleton random number generator.
0015  */
0016 class RNG {
0017
0018 private:
0019   std::mt19937 rng;
0020   std::uniform_real_distribution<> uniform_dist;
0021
0022   VECGEOM_FORCE_INLINE
0023   Precision GetUniform() { return uniform_dist(rng); }
0024
0025
0026 public:
0027   RNG() : rng(0), uniform_dist(0, 1) {}
0028
0029 public:
0030   void seed(unsigned long seed_val) { rng.seed(seed_val); }
0031
0032   /**
0033    * Access singleton instance.
0034    */
0035   static RNG &Instance()
0036   {
0037     static RNG instance;
0038     return instance;
0039   }
0040
0041   /**
0042    * @return Uniformly distributed floating point number between 0 and 1 unless
0043    *         range arguments are passed.
0044    */
0045   Precision uniform(const Precision min = 0., const Precision max = 1.) { return min + (max - min) * GetUniform(); }
0046
0047   int Poisson(const Precision lambda)
0048   {
0049     int k                  = 0;
0050     const Precision target = exp(-lambda);
0051     Precision p            = GetUniform();
0052     while (p < target) {
0053       p *= GetUniform();
0054       ++k;
0055     }
0056     return k;
0057   }
0058
0059   // interface for ROOT compatibility
0060   Precision Gaus(Precision ave = 0.0, Precision sig = 1.0) { return Gauss(ave, sig); }
0061
0062   Precision Gauss(Precision ave = 0.0, Precision sig = 1.0)
0063   {
0064     Precision x1, x2, w;
0065
0066     do {
0067       x1 = 2.0 * GetUniform() - 1.0;
0068       x2 = 2.0 * GetUniform() - 1.0;
0069       w  = x1 * x1 + x2 * x2;
0070     } while (w >= 1.0);
0071
0072     w = std::sqrt((-2.0 * std::log(w)) / w);
0073     return ave + (x1 * w * sig);
0074   }
0075
0076   /**
0077    * Uniformly distributed array of floating point number between 0 and 1 unless
0078    *         range arguments are passed.
0079    */
0080   void uniform_array(size_t n, Precision *array, const Precision min = 0., const Precision max = 1.)
0081   {
0082     for (size_t i = 0; i < n; ++i) {
0083       array[i] = min + (max - min) * GetUniform();
0084     }
0085   }
0086
0087 private:
0088   RNG(RNG const &);
0089   RNG &operator=(RNG const &);
0090 };
0091
0092 } // End global namespace
0093
0094 #endif // VECGEOM_BASE_RNG_H_