Geant4 Cross Reference |
1 /* 2 Code adapted from: 3 http://www.jstatsoft.org/v05/i08/ 4 5 6 Original disclaimer: 7 8 The ziggurat method for RNOR and REXP 9 Combine the code below with the main program in which you want 10 normal or exponential variates. Then use of RNOR in any expression 11 will provide a standard normal variate with mean zero, variance 1, 12 while use of REXP in any expression will provide an exponential variate 13 with density exp(-x),x>0. 14 Before using RNOR or REXP in your main, insert a command such as 15 zigset(86947731 ); 16 with your own choice of seed value>0, rather than 86947731. 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 Tsang, "The ziggurat method 19 for generating random variables", Journ. Statistical Software. 20 21 */ 22 23 #ifndef RandExpZiggurat_h 24 #define RandExpZiggurat_h 1 25 26 #include "CLHEP/Random/Random.h" 27 #include "CLHEP/Utility/memory.h" 28 #include "CLHEP/Utility/thread_local.h" 29 30 namespace CLHEP { 31 32 /** 33 * @author ATLAS 34 * @ingroup random 35 */ 36 class RandExpZiggurat : public HepRandom { 37 38 public: 39 40 inline RandExpZiggurat ( HepRandomEngine& anEngine, double mean=1.0 ); 41 inline RandExpZiggurat ( HepRandomEngine* anEngine, double mean=1.0 ); 42 // These constructors should be used to instantiate a RandExpZiggurat 43 // distribution object defining a local engine for it. 44 // The static generator will be skipped using the non-static methods 45 // defined below. 46 // If the engine is passed by pointer the corresponding engine object 47 // will be deleted by the RandExpZiggurat destructor. 48 // If the engine is passed by reference the corresponding engine object 49 // will not be deleted by the RandExpZiggurat destructor. 50 51 virtual ~RandExpZiggurat(); 52 // Destructor 53 54 // Static methods to shoot random values using the static generator 55 56 static float shoot() {return shoot(HepRandom::getTheEngine());} 57 static float shoot( float mean ) {return shoot(HepRandom::getTheEngine(),mean);} 58 59 /* ENGINE IS INTRINSIC FLOAT 60 static double shoot() {return shoot(HepRandom::getTheEngine());} 61 static double shoot( double mean ) {return shoot(HepRandom::getTheEngine(),mean);} 62 */ 63 64 static void shootArray ( const int size, float* vect, float mean=1.0 ); 65 static void shootArray ( const int size, double* vect, double mean=1.0 ); 66 67 // Static methods to shoot random values using a given engine 68 // by-passing the static generator. 69 70 static inline float shoot( HepRandomEngine* anEngine ) {return ziggurat_REXP(anEngine);} 71 static inline float shoot( HepRandomEngine* anEngine, float mean ) {return shoot(anEngine)*mean;} 72 73 /* ENGINE IS INTRINSIC FLOAT 74 static inline double shoot( HepRandomEngine* anEngine ) {return ziggurat_REXP(anEngine);} 75 76 static inline double shoot( HepRandomEngine* anEngine, double mean ) {return shoot(anEngine)*mean;} 77 */ 78 79 static void shootArray ( HepRandomEngine* anEngine, const int size, float* vect, float mean=1.0 ); 80 static void shootArray ( HepRandomEngine* anEngine, const int size, double* vect, double mean=1.0 ); 81 82 // Methods using the localEngine to shoot random values, by-passing 83 // the static generator. 84 85 inline float fire() {return fire(float(defaultMean));} 86 inline float fire( float mean ) {return ziggurat_REXP(localEngine.get())*mean;} 87 88 /* ENGINE IS INTRINSIC FLOAT 89 inline double fire() {return fire(defaultMean);} 90 inline double fire( double mean ) {return ziggurat_REXP(localEngine.get())*mean;} 91 */ 92 93 void fireArray ( const int size, float* vect ); 94 void fireArray ( const int size, double* vect ); 95 void fireArray ( const int size, float* vect, float mean ); 96 void fireArray ( const int size, double* vect, double mean ); 97 98 virtual double operator()(); 99 inline float operator()( float mean ) {return fire( mean );} 100 101 // Save and restore to/from streams 102 103 std::ostream & put ( std::ostream & os ) const; 104 std::istream & get ( std::istream & is ); 105 106 std::string name() const; 107 HepRandomEngine & engine(); 108 109 static std::string distributionName() {return "RandExpZiggurat";} 110 // Provides the name of this distribution class 111 112 static bool ziggurat_init(); 113 protected: 114 ////////////////////////// 115 // Ziggurat Original code: 116 ////////////////////////// 117 //static unsigned long jz,jsr=123456789; 118 // 119 //#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr) 120 //#define UNI (.5 + (signed) SHR3*.2328306e-9) 121 //#define IUNI SHR3 122 // 123 //static long hz; 124 //static unsigned long iz, kn[128], ke[256]; 125 //static float wn[128],fn[128], we[256],fe[256]; 126 // 127 //#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz)<kn[iz])? hz*wn[iz] : nfix()) 128 //#define REXP (jz=SHR3, iz=jz&255, ( jz <ke[iz])? jz*we[iz] : efix()) 129 130 static CLHEP_THREAD_LOCAL unsigned long kn[128], ke[256]; 131 static CLHEP_THREAD_LOCAL float wn[128],fn[128], we[256],fe[256]; 132 133 static CLHEP_THREAD_LOCAL bool ziggurat_is_init; 134 135 static inline unsigned long ziggurat_SHR3(HepRandomEngine* anEngine) {return (unsigned int)(*anEngine);} 136 static inline float ziggurat_UNI(HepRandomEngine* anEngine) {return float(anEngine->flat());} 137 static inline float ziggurat_REXP(HepRandomEngine* anEngine) { 138 if(!ziggurat_is_init) ziggurat_init(); 139 unsigned long jz=ziggurat_SHR3(anEngine); 140 unsigned long iz=jz&255; 141 return (jz<ke[iz]) ? jz*we[iz] : ziggurat_efix(jz,anEngine); 142 } 143 static float ziggurat_efix(unsigned long jz,HepRandomEngine* anEngine); 144 145 private: 146 147 // Private copy constructor. Defining it here disallows use. 148 RandExpZiggurat(const RandExpZiggurat& d); 149 150 std::shared_ptr<HepRandomEngine> localEngine; 151 double defaultMean; 152 }; 153 154 } // namespace CLHEP 155 156 namespace CLHEP { 157 158 inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine & anEngine, double mean ) : localEngine(&anEngine, do_nothing_deleter()), defaultMean(mean) 159 { 160 } 161 162 inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine * anEngine, double mean ) : localEngine(anEngine), defaultMean(mean) 163 { 164 } 165 166 } // namespace CLHEP 167 168 #endif 169