include/G4HepEm/G4HepEmRandomEngine.hh

0001
0002 #ifndef G4HepEmRandomEngine_HH
0003 #define G4HepEmRandomEngine_HH
0004
0005 #include "G4HepEmMacros.hh"
0006 #include "G4HepEmMath.hh"
0007 #include "G4HepEmConstants.hh"
0008
0009 #include <cmath>
0010
0011 /**
0012  * @file    G4HepEmRandomEngine.hh
0013  * @class   G4HepEmRandomEngine
0014  * @author  J. Hahnfeld
0015  * @date    2021
0016  *
0017  * A simple abstraction for a random number engine.
0018  *
0019  * Holds a reference to the real engine and two member functions to use the
0020  * engine to generate one random number or fill an array with a given size, respectively.
0021  *
0022  * When G4HepEm is compiled with support for Geant4, the reference must be a pointer
0023  * to an instance of `CLHEP::HepRandomEngine` for host-side use. For device-side use,
0024  * the user must implement a suitable reference _and_ compile/link `__device__` implementations
0025  * for the `flat` and `flatArray` member functions.
0026  *
0027  * For G4HepEm built in standalone mode without Geant4 support, the user must compile and
0028  * link in both host- and device- side implementations for the engine and member functions.
0029  */
0030 class G4HepEmRandomEngine final {
0031 public:
0032   G4HepEmHostDevice
0033   G4HepEmRandomEngine(void *object)
0034     : fObject(object), fIsGauss(false), fGauss(0.) { }
0035
0036   G4HepEmHostDevice
0037   void SetObject(void* p) { fObject = p; }
0038
0039   /** Return a random number uniformly distributed between 0 and 1.
0040    */
0041   G4HepEmHostDevice
0042   double flat();
0043   /** Fill elements of array with random numbers uniformly distributed between 0 and 1.
0044    *
0045    *  @param [in] size Number of elements in `vect` input array
0046    *  @param [in][out] vect Array to fill with random numbers
0047    *  @pre `size` must be less than or equal to the number of elements in `vect`
0048    */
0049   G4HepEmHostDevice
0050   void flatArray(const int size, double* vect);
0051
0052   G4HepEmHostDevice
0053   double Gauss(const double mean, const double stDev) {
0054     if (fIsGauss) {
0055       fIsGauss = false;
0056       return fGauss*stDev+mean;
0057     }
0058     double rnd[2];
0059     double r, v1, v2;
0060     do {
0061       flatArray(2, rnd);
0062       v1 = 2.*rnd[0] - 1.;
0063       v2 = 2.*rnd[1] - 1.;
0064       r = v1*v1 + v2*v2;
0065     } while ( r > 1.);
0066     const double fac = std::sqrt(-2.*G4HepEmLog(r)/r);
0067     fGauss   = v1*fac;
0068     fIsGauss = true;
0069     return v2*fac*stDev+mean;
0070   }
0071
0072   G4HepEmHostDevice
0073   void DiscardGauss() { fIsGauss = false; }
0074
0075
0076   G4HepEmHostDevice
0077   int Poisson(double mean) {
0078     const int   border = 16;
0079     const double limit = 2.E+9;
0080
0081     int number = 0;
0082     if(mean <= border) {
0083       const double position = flat();
0084       double poissonValue   = G4HepEmExp(-mean);
0085       double poissonSum     = poissonValue;
0086       while(poissonSum <= position) {
0087         ++number;
0088         poissonValue *= mean/number;
0089         poissonSum   += poissonValue;
0090       }
0091       return number;
0092     }  // the case of mean <= 16
0093     //
0094     double rnd[2];
0095     flatArray(2, rnd);
0096     const double t = std::sqrt(-2.*G4HepEmLog(rnd[0])) * std::cos(k2Pi*rnd[1]);
0097     double value = mean + t*std::sqrt(mean) + 0.5;
0098     return value < 0.     ?  0 :
0099            value >= limit ? static_cast<int>(limit) : static_cast<int>(value);
0100   }
0101
0102
0103 private:
0104   void *fObject;
0105
0106   bool fIsGauss;
0107   double fGauss;
0108 };
0109
0110 #endif // G4HepEmRandomEngine_HH