| 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 #ifndef RandExpZiggurat_h 23 #ifndef RandExpZiggurat_h
24 #define RandExpZiggurat_h 1 24 #define RandExpZiggurat_h 1
25 25
26 #include "CLHEP/Random/Random.h" 26 #include "CLHEP/Random/Random.h"
27 #include "CLHEP/Utility/memory.h" <<
28 #include "CLHEP/Utility/thread_local.h" 27 #include "CLHEP/Utility/thread_local.h"
29 28
30 namespace CLHEP { 29 namespace CLHEP {
31 30
32 /** 31 /**
33 * @author ATLAS 32 * @author ATLAS
34 * @ingroup random 33 * @ingroup random
35 */ 34 */
36 class RandExpZiggurat : public HepRandom { 35 class RandExpZiggurat : public HepRandom {
37 36
38 public: 37 public:
39 38
40 inline RandExpZiggurat ( HepRandomEngine& an 39 inline RandExpZiggurat ( HepRandomEngine& anEngine, double mean=1.0 );
41 inline RandExpZiggurat ( HepRandomEngine* an 40 inline RandExpZiggurat ( HepRandomEngine* anEngine, double mean=1.0 );
42 // These constructors should be used to inst 41 // These constructors should be used to instantiate a RandExpZiggurat
43 // distribution object defining a local engi 42 // distribution object defining a local engine for it.
44 // The static generator will be skipped usin 43 // The static generator will be skipped using the non-static methods
45 // defined below. 44 // defined below.
46 // If the engine is passed by pointer the co 45 // If the engine is passed by pointer the corresponding engine object
47 // will be deleted by the RandExpZiggurat de 46 // will be deleted by the RandExpZiggurat destructor.
48 // If the engine is passed by reference the 47 // If the engine is passed by reference the corresponding engine object
49 // will not be deleted by the RandExpZiggura 48 // will not be deleted by the RandExpZiggurat destructor.
50 49
51 virtual ~RandExpZiggurat(); 50 virtual ~RandExpZiggurat();
52 // Destructor 51 // Destructor
53 52
54 // Static methods to shoot random values usi 53 // Static methods to shoot random values using the static generator
55 54
56 static float shoot() {return shoot(HepRandom << 55 static float shoot() {return shoot(HepRandom::getTheEngine());};
57 static float shoot( float mean ) {return sho << 56 static float shoot( float mean ) {return shoot(HepRandom::getTheEngine(),mean);};
58 57
59 /* ENGINE IS INTRINSIC FLOAT 58 /* ENGINE IS INTRINSIC FLOAT
60 static double shoot() {return shoot(HepRando << 59 static double shoot() {return shoot(HepRandom::getTheEngine());};
61 static double shoot( double mean ) {return s << 60 static double shoot( double mean ) {return shoot(HepRandom::getTheEngine(),mean);};
62 */ 61 */
63 62
64 static void shootArray ( const int size, flo 63 static void shootArray ( const int size, float* vect, float mean=1.0 );
65 static void shootArray ( const int size, dou 64 static void shootArray ( const int size, double* vect, double mean=1.0 );
66 65
67 // Static methods to shoot random values us 66 // Static methods to shoot random values using a given engine
68 // by-passing the static generator. 67 // by-passing the static generator.
69 68
70 static inline float shoot( HepRandomEngine* << 69 static inline float shoot( HepRandomEngine* anEngine ) {return ziggurat_REXP(anEngine);};
71 static inline float shoot( HepRandomEngine* << 70 static inline float shoot( HepRandomEngine* anEngine, float mean ) {return shoot(anEngine)*mean;};
72 71
73 /* ENGINE IS INTRINSIC FLOAT 72 /* ENGINE IS INTRINSIC FLOAT
74 static inline double shoot( HepRandomEngine* << 73 static inline double shoot( HepRandomEngine* anEngine ) {return ziggurat_REXP(anEngine);};
75 74
76 static inline double shoot( HepRandomEngine* << 75 static inline double shoot( HepRandomEngine* anEngine, double mean ) {return shoot(anEngine)*mean;};
77 */ 76 */
78 77
79 static void shootArray ( HepRandomEngine* an 78 static void shootArray ( HepRandomEngine* anEngine, const int size, float* vect, float mean=1.0 );
80 static void shootArray ( HepRandomEngine* an 79 static void shootArray ( HepRandomEngine* anEngine, const int size, double* vect, double mean=1.0 );
81 80
82 // Methods using the localEngine to shoot r 81 // Methods using the localEngine to shoot random values, by-passing
83 // the static generator. 82 // the static generator.
84 83
85 inline float fire() {return fire(float(defau << 84 inline float fire() {return fire(defaultMean);};
86 inline float fire( float mean ) {return zigg << 85 inline float fire( float mean ) {return ziggurat_REXP(localEngine)*mean;};
87 86
88 /* ENGINE IS INTRINSIC FLOAT 87 /* ENGINE IS INTRINSIC FLOAT
89 inline double fire() {return fire(defaultMea << 88 inline double fire() {return fire(defaultMean);};
90 inline double fire( double mean ) {return zi << 89 inline double fire( double mean ) {return ziggurat_REXP(localEngine)*mean;};
91 */ 90 */
92 91
93 void fireArray ( const int size, float* vect 92 void fireArray ( const int size, float* vect );
94 void fireArray ( const int size, double* vec 93 void fireArray ( const int size, double* vect );
95 void fireArray ( const int size, float* vect 94 void fireArray ( const int size, float* vect, float mean );
96 void fireArray ( const int size, double* vec 95 void fireArray ( const int size, double* vect, double mean );
97 96
98 virtual double operator()(); 97 virtual double operator()();
99 inline float operator()( float mean ) {retur << 98 inline float operator()( float mean ) {return fire( mean );};
100 99
101 // Save and restore to/from streams 100 // Save and restore to/from streams
102 101
103 std::ostream & put ( std::ostream & os ) con 102 std::ostream & put ( std::ostream & os ) const;
104 std::istream & get ( std::istream & is ); 103 std::istream & get ( std::istream & is );
105 104
106 std::string name() const; 105 std::string name() const;
107 HepRandomEngine & engine(); 106 HepRandomEngine & engine();
108 107
109 static std::string distributionName() {retur 108 static std::string distributionName() {return "RandExpZiggurat";}
110 // Provides the name of this distribution cl 109 // Provides the name of this distribution class
111 110
112 static bool ziggurat_init(); 111 static bool ziggurat_init();
113 protected: 112 protected:
114 ////////////////////////// 113 //////////////////////////
115 // Ziggurat Original code: 114 // Ziggurat Original code:
116 ////////////////////////// 115 //////////////////////////
117 //static unsigned long jz,jsr=123456789; 116 //static unsigned long jz,jsr=123456789;
118 // 117 //
119 //#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^ 118 //#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
120 //#define UNI (.5 + (signed) SHR3*.2328306e- 119 //#define UNI (.5 + (signed) SHR3*.2328306e-9)
121 //#define IUNI SHR3 120 //#define IUNI SHR3
122 // 121 //
123 //static long hz; 122 //static long hz;
124 //static unsigned long iz, kn[128], ke[256]; 123 //static unsigned long iz, kn[128], ke[256];
125 //static float wn[128],fn[128], we[256],fe[2 124 //static float wn[128],fn[128], we[256],fe[256];
126 // 125 //
127 //#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz 126 //#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz)<kn[iz])? hz*wn[iz] : nfix())
128 //#define REXP (jz=SHR3, iz=jz&255, ( jz 127 //#define REXP (jz=SHR3, iz=jz&255, ( jz <ke[iz])? jz*we[iz] : efix())
129 128
130 static CLHEP_THREAD_LOCAL unsigned long kn[1 129 static CLHEP_THREAD_LOCAL unsigned long kn[128], ke[256];
131 static CLHEP_THREAD_LOCAL float wn[128],fn[1 130 static CLHEP_THREAD_LOCAL float wn[128],fn[128], we[256],fe[256];
132 131
133 static CLHEP_THREAD_LOCAL bool ziggurat_is_i 132 static CLHEP_THREAD_LOCAL bool ziggurat_is_init;
134 133
135 static inline unsigned long ziggurat_SHR3(He << 134 static inline unsigned long ziggurat_SHR3(HepRandomEngine* anEngine) {return (unsigned int)(*anEngine);};
136 static inline float ziggurat_UNI(HepRandomEn << 135 static inline float ziggurat_UNI(HepRandomEngine* anEngine) {return anEngine->flat();};
137 static inline float ziggurat_REXP(HepRandomE 136 static inline float ziggurat_REXP(HepRandomEngine* anEngine) {
138 if(!ziggurat_is_init) ziggurat_init(); 137 if(!ziggurat_is_init) ziggurat_init();
139 unsigned long jz=ziggurat_SHR3(anEngine); 138 unsigned long jz=ziggurat_SHR3(anEngine);
140 unsigned long iz=jz&255; 139 unsigned long iz=jz&255;
141 return (jz<ke[iz]) ? jz*we[iz] : ziggurat_ 140 return (jz<ke[iz]) ? jz*we[iz] : ziggurat_efix(jz,anEngine);
142 } << 141 };
143 static float ziggurat_efix(unsigned long jz, 142 static float ziggurat_efix(unsigned long jz,HepRandomEngine* anEngine);
144 143
145 private: 144 private:
146 145
147 // Private copy constructor. Defining it her 146 // Private copy constructor. Defining it here disallows use.
148 RandExpZiggurat(const RandExpZiggurat& d); 147 RandExpZiggurat(const RandExpZiggurat& d);
149 148
150 std::shared_ptr<HepRandomEngine> localEngine << 149 HepRandomEngine* localEngine;
>> 150 bool deleteEngine;
151 double defaultMean; 151 double defaultMean;
152 }; 152 };
153 153
154 } // namespace CLHEP 154 } // namespace CLHEP
155 155
156 namespace CLHEP { 156 namespace CLHEP {
157 157
158 inline RandExpZiggurat::RandExpZiggurat(HepRan << 158 inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine & anEngine, double mean ) : localEngine(&anEngine), deleteEngine(false), defaultMean(mean)
159 { 159 {
160 } 160 }
161 161
162 inline RandExpZiggurat::RandExpZiggurat(HepRan << 162 inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine * anEngine, double mean ) : localEngine(anEngine), deleteEngine(true), defaultMean(mean)
163 { 163 {
164 } 164 }
165 165
166 } // namespace CLHEP 166 } // namespace CLHEP
167 167
168 #endif 168 #endif
169 169