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