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 #ifdef VECCORE_CUDA
0010 #include <cuda.h>
0011 #include <curand_kernel.h>
0012 #else
0013 #include <random>
0014 #endif
0015
0016 namespace vecgeom {
0017 inline namespace VECGEOM_IMPL_NAMESPACE {
0018
0019 #ifdef VECCORE_CUDA
0020 class RNG;
0021
0022 // Emulating static class member ..
0023 namespace RNGvar {
0024 extern VECCORE_ATT_DEVICE unsigned long gMaxInstance;
0025 extern VECCORE_ATT_DEVICE RNG **gInstances;
0026 }
0027 #endif
0028
0029 /**
0030  * @brief Singleton random number generator.
0031  */
0032 class RNG {
0033
0034 private:
0035 #ifdef VECCORE_CUDA
0036
0037 #ifdef __CUDA_ARCH__
0038   curandState fState;
0039 #else
0040 // Using rand in C++03
0041 #endif
0042
0043   VECCORE_ATT_HOST_DEVICE
0044   VECGEOM_FORCE_INLINE
0045   Precision GetUniform()
0046   {
0047 #ifdef __CUDA_ARCH__
0048     return curand_uniform(&fState);
0049 #else
0050     return (Precision)rand() / RAND_MAX;
0051 #endif
0052   }
0053
0054 #else
0055
0056   std::mt19937 rng;
0057   std::uniform_real_distribution<> uniform_dist;
0058
0059   VECGEOM_FORCE_INLINE
0060   Precision GetUniform() { return uniform_dist(rng); }
0061
0062 #endif
0063
0064 public:
0065 #ifdef VECCORE_CUDA
0066   // The state should really be 'thread' specific
0067   VECCORE_ATT_HOST_DEVICE
0068   RNG()
0069   {
0070 #ifdef __CUDA_ARCH__
0071     curand_init(0 /*seed */, 0 /* subsequence */, 0 /* offset */, &fState);
0072 #else
0073 // using rand in C++03
0074 #endif
0075   }
0076 #else
0077   RNG() : rng(0), uniform_dist(0, 1) {}
0078 #endif
0079
0080 public:
0081 /**
0082  * Init thread specific singleton instance.
0083  */
0084 #ifdef __CUDA_ARCH__
0085   VECCORE_ATT_DEVICE
0086   static void InitInstances(unsigned long nthreads)
0087   {
0088     unsigned int tid = (threadIdx.x + blockIdx.x * blockDim.x);
0089
0090     if (tid == 0) {
0091       RNGvar::gMaxInstance = nthreads;
0092       RNGvar::gInstances   = new RNG *[nthreads];
0093     }
0094     __syncthreads();
0095
0096     for (int i = tid; i < nthreads; i += blockDim.x * gridDim.x) {
0097       RNGvar::gInstances[i] = new RNG;
0098     }
0099   }
0100 #endif
0101
0102 #ifndef VECCORE_CUDA
0103   void seed(unsigned long seed_val) { rng.seed(seed_val); }
0104 #endif
0105   /**
0106    * Access singleton instance.
0107    */
0108   VECCORE_ATT_HOST_DEVICE
0109   static RNG &Instance()
0110   {
0111 #ifdef __CUDA_ARCH__
0112     unsigned int tid = (threadIdx.x + blockIdx.x * blockDim.x);
0113     if (tid < RNGvar::gMaxInstance)
0114       return *(RNGvar::gInstances[tid]);
0115     else
0116       return *(new RNG);
0117 #else
0118     static RNG instance;
0119     return instance;
0120 #endif
0121   }
0122
0123   /**
0124    * @return Uniformly distributed floating point number between 0 and 1 unless
0125    *         range arguments are passed.
0126    */
0127   VECCORE_ATT_HOST_DEVICE
0128   VECGEOM_FORCE_INLINE
0129   Precision uniform(const Precision min = 0., const Precision max = 1.) { return min + (max - min) * GetUniform(); }
0130
0131   VECCORE_ATT_HOST_DEVICE
0132   VECGEOM_FORCE_INLINE
0133   int Poisson(const Precision lambda)
0134   {
0135     int k                  = 0;
0136     const Precision target = exp(-lambda);
0137     Precision p            = GetUniform();
0138     while (p < target) {
0139       p *= GetUniform();
0140       ++k;
0141     }
0142     return k;
0143   }
0144
0145   // interface for ROOT compatibility
0146   VECCORE_ATT_HOST_DEVICE
0147   VECGEOM_FORCE_INLINE
0148   Precision Gaus(Precision ave = 0.0, Precision sig = 1.0) { return Gauss(ave, sig); }
0149
0150   VECCORE_ATT_HOST_DEVICE
0151   VECGEOM_FORCE_INLINE
0152   Precision Gauss(Precision ave = 0.0, Precision sig = 1.0)
0153   {
0154     Precision x1, x2, w;
0155
0156     do {
0157       x1 = 2.0 * GetUniform() - 1.0;
0158       x2 = 2.0 * GetUniform() - 1.0;
0159       w  = x1 * x1 + x2 * x2;
0160     } while (w >= 1.0);
0161
0162     w = std::sqrt((-2.0 * std::log(w)) / w);
0163     return ave + (x1 * w * sig);
0164   }
0165
0166   /**
0167    * Uniformly distributed array of floating point number between 0 and 1 unless
0168    *         range arguments are passed.
0169    */
0170   VECCORE_ATT_HOST_DEVICE
0171   VECGEOM_FORCE_INLINE
0172   void uniform_array(size_t n, Precision *array, const Precision min = 0., const Precision max = 1.)
0173   {
0174     for (size_t i = 0; i < n; ++i) {
0175       array[i] = min + (max - min) * GetUniform();
0176     }
0177   }
0178
0179 private:
0180   RNG(RNG const &);
0181   RNG &operator=(RNG const &);
0182 };
0183 }
0184 } // End global namespace
0185
0186 #endif // VECGEOM_BASE_RNG_H_