Geant4 Cross Reference |
1 /* 1 /* 2 Code adapted from: 2 Code adapted from: 3 http://www.jstatsoft.org/v05/i08/ 3 http://www.jstatsoft.org/v05/i08/ 4 4 5 5 6 Original disclaimer: 6 Original disclaimer: 7 7 8 The ziggurat method for RNOR and REXP 8 The ziggurat method for RNOR and REXP 9 Combine the code below with the main program i 9 Combine the code below with the main program in which you want 10 normal or exponential variates. Then use of 10 normal or exponential variates. Then use of RNOR in any expression 11 will provide a standard normal variate with me 11 will provide a standard normal variate with mean zero, variance 1, 12 while use of REXP in any expression will provi 12 while use of REXP in any expression will provide an exponential variate 13 with density exp(-x),x>0. 13 with density exp(-x),x>0. 14 Before using RNOR or REXP in your main, insert 14 Before using RNOR or REXP in your main, insert a command such as 15 zigset(86947731 ); 15 zigset(86947731 ); 16 with your own choice of seed value>0, rather t 16 with your own choice of seed value>0, rather than 86947731. 17 (If you do not invoke zigset(...) you will get 17 (If you do not invoke zigset(...) you will get all zeros for RNOR and REXP.) 18 For details of the method, see Marsaglia and T 18 For details of the method, see Marsaglia and Tsang, "The ziggurat method 19 for generating random variables", Journ. Stati 19 for generating random variables", Journ. Statistical Software. 20 20 21 */ 21 */ 22 22 23 23 24 #ifndef RandGaussZiggurat_h 24 #ifndef RandGaussZiggurat_h 25 #define RandGaussZiggurat_h 1 25 #define RandGaussZiggurat_h 1 26 26 27 #include "cmath" 27 #include "cmath" 28 #include "CLHEP/Random/RandGauss.h" 28 #include "CLHEP/Random/RandGauss.h" 29 #include "CLHEP/Utility/thread_local.h" 29 #include "CLHEP/Utility/thread_local.h" 30 30 31 namespace CLHEP { 31 namespace CLHEP { 32 32 33 /** 33 /** 34 * @author ATLAS 34 * @author ATLAS 35 * @ingroup random 35 * @ingroup random 36 */ 36 */ 37 class RandGaussZiggurat : public RandGauss { 37 class RandGaussZiggurat : public RandGauss { 38 38 39 public: 39 public: 40 40 41 inline RandGaussZiggurat ( HepRandomEngine& 41 inline RandGaussZiggurat ( HepRandomEngine& anEngine, double mean=0.0, double stdDev=1.0 ); 42 inline RandGaussZiggurat ( HepRandomEngine* 42 inline RandGaussZiggurat ( HepRandomEngine* anEngine, double mean=0.0, double stdDev=1.0 ); 43 43 44 // Destructor 44 // Destructor 45 virtual ~RandGaussZiggurat(); 45 virtual ~RandGaussZiggurat(); 46 46 47 // Static methods to shoot random values usi 47 // Static methods to shoot random values using the static generator 48 48 49 static inline float shoot() {return ziggurat << 49 static inline float shoot() {return ziggurat_RNOR(HepRandom::getTheEngine());}; 50 static inline float shoot( float mean, float << 50 static inline float shoot( float mean, float stdDev ) {return shoot()*stdDev + mean;}; 51 51 52 static void shootArray ( const int size, flo 52 static void shootArray ( const int size, float* vect, float mean=0.0, float stdDev=1.0 ); 53 static void shootArray ( const int size, dou 53 static void shootArray ( const int size, double* vect, double mean=0.0, double stdDev=1.0 ); 54 54 55 // Static methods to shoot random values us 55 // Static methods to shoot random values using a given engine 56 // by-passing the static generator. 56 // by-passing the static generator. 57 57 58 static inline float shoot( HepRandomEngine* << 58 static inline float shoot( HepRandomEngine* anotherEngine ) {return ziggurat_RNOR(anotherEngine);}; 59 static inline float shoot( HepRandomEngine* << 59 static inline float shoot( HepRandomEngine* anotherEngine, float mean, float stdDev ) {return shoot(anotherEngine)*stdDev + mean;}; 60 60 61 static void shootArray ( HepRandomEngine* an 61 static void shootArray ( HepRandomEngine* anotherEngine, const int size, float* vect, float mean=0.0, float stdDev=1.0 ); 62 static void shootArray ( HepRandomEngine* an 62 static void shootArray ( HepRandomEngine* anotherEngine, const int size, double* vect, double mean=0.0, double stdDev=1.0 ); 63 63 64 // Instance methods using the localEngine t 64 // Instance methods using the localEngine to instead of the static 65 // generator, and the default mean and stdD 65 // generator, and the default mean and stdDev established at construction 66 66 67 inline float fire() {return float(ziggurat_R << 67 inline float fire() {return ziggurat_RNOR(localEngine.get()) * defaultStdDev + defaultMean;}; 68 68 69 inline float fire ( float mean, float stdDev << 69 inline float fire ( float mean, float stdDev ) {return ziggurat_RNOR(localEngine.get()) * stdDev + mean;}; 70 70 71 void fireArray ( const int size, float* vec 71 void fireArray ( const int size, float* vect); 72 void fireArray ( const int size, double* ve 72 void fireArray ( const int size, double* vect); 73 void fireArray ( const int size, float* vec 73 void fireArray ( const int size, float* vect, float mean, float stdDev ); 74 void fireArray ( const int size, double* ve 74 void fireArray ( const int size, double* vect, double mean, double stdDev ); 75 75 76 virtual double operator()(); 76 virtual double operator()(); 77 virtual double operator()( double mean, doub 77 virtual double operator()( double mean, double stdDev ); 78 78 79 // Save and restore to/from streams 79 // Save and restore to/from streams 80 80 81 std::ostream & put ( std::ostream & os ) con 81 std::ostream & put ( std::ostream & os ) const; 82 std::istream & get ( std::istream & is ); 82 std::istream & get ( std::istream & is ); 83 83 84 std::string name() const; 84 std::string name() const; 85 HepRandomEngine & engine(); 85 HepRandomEngine & engine(); 86 86 87 static std::string distributionName() {retur 87 static std::string distributionName() {return "RandGaussZiggurat";} 88 // Provides the name of this distribution cl 88 // Provides the name of this distribution class 89 89 90 static bool ziggurat_init(); 90 static bool ziggurat_init(); 91 protected: 91 protected: 92 ////////////////////////// 92 ////////////////////////// 93 // Ziggurat Original code: 93 // Ziggurat Original code: 94 ////////////////////////// 94 ////////////////////////// 95 //static unsigned long jz,jsr=123456789; 95 //static unsigned long jz,jsr=123456789; 96 // 96 // 97 //#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^ 97 //#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr) 98 //#define UNI (.5 + (signed) SHR3*.2328306e- 98 //#define UNI (.5 + (signed) SHR3*.2328306e-9) 99 //#define IUNI SHR3 99 //#define IUNI SHR3 100 // 100 // 101 //static long hz; 101 //static long hz; 102 //static unsigned long iz, kn[128], ke[256]; 102 //static unsigned long iz, kn[128], ke[256]; 103 //static float wn[128],fn[128], we[256],fe[2 103 //static float wn[128],fn[128], we[256],fe[256]; 104 // 104 // 105 //#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz 105 //#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz)<kn[iz])? hz*wn[iz] : nfix()) 106 //#define REXP (jz=SHR3, iz=jz&255, ( jz 106 //#define REXP (jz=SHR3, iz=jz&255, ( jz <ke[iz])? jz*we[iz] : efix()) 107 107 108 static CLHEP_THREAD_LOCAL unsigned long kn[1 108 static CLHEP_THREAD_LOCAL unsigned long kn[128], ke[256]; 109 static CLHEP_THREAD_LOCAL float wn[128],fn[1 109 static CLHEP_THREAD_LOCAL float wn[128],fn[128], we[256],fe[256]; 110 110 111 static CLHEP_THREAD_LOCAL bool ziggurat_is_i 111 static CLHEP_THREAD_LOCAL bool ziggurat_is_init; 112 112 113 static inline unsigned long ziggurat_SHR3(He << 113 static inline unsigned long ziggurat_SHR3(HepRandomEngine* anEngine) {return (unsigned int)(*anEngine);}; 114 static inline float ziggurat_UNI(HepRandomEn << 114 static inline float ziggurat_UNI(HepRandomEngine* anEngine) {return anEngine->flat();}; 115 static inline float ziggurat_RNOR(HepRandomE 115 static inline float ziggurat_RNOR(HepRandomEngine* anEngine) { 116 if(!ziggurat_is_init) ziggurat_init(); 116 if(!ziggurat_is_init) ziggurat_init(); 117 long hz=(signed)ziggurat_SHR3(anEngine); 117 long hz=(signed)ziggurat_SHR3(anEngine); 118 unsigned long iz=hz&127; 118 unsigned long iz=hz&127; 119 return ((unsigned long)std::abs(hz)<kn[iz] 119 return ((unsigned long)std::abs(hz)<kn[iz]) ? hz*wn[iz] : ziggurat_nfix(hz,anEngine); 120 } << 120 }; 121 static float ziggurat_nfix(long hz,HepRandom 121 static float ziggurat_nfix(long hz,HepRandomEngine* anEngine); 122 122 123 private: 123 private: 124 124 125 // Private copy constructor. Defining it her 125 // Private copy constructor. Defining it here disallows use. 126 RandGaussZiggurat(const RandGaussZiggurat& d 126 RandGaussZiggurat(const RandGaussZiggurat& d); 127 127 128 // All the engine info, and the default mean 128 // All the engine info, and the default mean and sigma, are in the RandGauss 129 // base class. 129 // base class. 130 130 131 }; 131 }; 132 132 133 } // namespace CLHEP 133 } // namespace CLHEP 134 134 135 namespace CLHEP { 135 namespace CLHEP { 136 136 137 RandGaussZiggurat::RandGaussZiggurat(HepRandom 137 RandGaussZiggurat::RandGaussZiggurat(HepRandomEngine & anEngine, double mean,double stdDev ): RandGauss(anEngine, mean, stdDev) 138 { 138 { 139 } 139 } 140 140 141 RandGaussZiggurat::RandGaussZiggurat(HepRandom 141 RandGaussZiggurat::RandGaussZiggurat(HepRandomEngine * anEngine, double mean,double stdDev ): RandGauss(anEngine, mean, stdDev) 142 { 142 { 143 } 143 } 144 144 145 } // namespace CLHEP 145 } // namespace CLHEP 146 146 147 #endif 147 #endif 148 148