ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RandGaussZiggurat.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RandGaussZiggurat.cc
3 #include <iostream>
4 #include <cmath> // for std::log()
5 
6 namespace CLHEP {
7 
11 
13 
15 }
16 
17 std::string RandGaussZiggurat::name() const
18 {
19  return "RandGaussZiggurat";
20 }
21 
23 {
24  const double rzm1 = 2147483648.0, rzm2 = 4294967296.;
25  double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
26  double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
27  int i;
28 
29 /* Set up tables for RNOR */
30  q=vn/std::exp(-.5*dn*dn);
31  kn[0]=(unsigned long)((dn/q)*rzm1);
32  kn[1]=0;
33 
34  wn[0]=q/rzm1;
35  wn[127]=dn/rzm1;
36 
37  fn[0]=1.;
38  fn[127]=std::exp(-.5*dn*dn);
39 
40  for(i=126;i>=1;i--) {
41  dn=std::sqrt(-2.*std::log(vn/dn+std::exp(-.5*dn*dn)));
42  kn[i+1]=(unsigned long)((dn/tn)*rzm1);
43  tn=dn;
44  fn[i]=std::exp(-.5*dn*dn);
45  wn[i]=dn/rzm1;
46  }
47 
48 /* Set up tables for REXP */
49  q = ve/std::exp(-de);
50  ke[0]=(unsigned long)((de/q)*rzm2);
51  ke[1]=0;
52 
53  we[0]=q/rzm2;
54  we[255]=de/rzm2;
55 
56  fe[0]=1.;
57  fe[255]=std::exp(-de);
58 
59  for(i=254;i>=1;i--) {
60  de=-std::log(ve/de+std::exp(-de));
61  ke[i+1]= (unsigned long)((de/te)*rzm2);
62  te=de;
63  fe[i]=std::exp(-de);
64  we[i]=de/rzm2;
65  }
66  ziggurat_is_init=true;
67 
68  //std::cout<<"Done RandGaussZiggurat::ziggurat_init()"<<std::endl;
69 
70  return true;
71 }
72 
74 {
76  const float r = 3.442620f; /* The start of the right tail */
77  float x, y;
78  unsigned long iz=hz&127;
79  for(;;)
80  {
81  x=hz*wn[iz]; /* iz==0, handles the base strip */
82  if(iz==0) {
83  do {
84  /* change to (1.0 - UNI) as argument to std::log(), because CLHEP generates [0,1),
85  while the original UNI generates (0,1] */
86  x=-std::log(1.0 - ziggurat_UNI(anEngine))*0.2904764; /* .2904764 is 1/r */
87  y=-std::log(1.0 - ziggurat_UNI(anEngine));
88  } while(y+y<x*x);
89  return (hz>0)? r+x : -r-x;
90  }
91  /* iz>0, handle the wedges of other strips */
92  if( fn[iz]+(1.0 - ziggurat_UNI(anEngine))*(fn[iz-1]-fn[iz]) < std::exp(-.5*x*x) ) return x;
93 
94  /* initiate, try to exit for(;;) for loop*/
95  hz=(signed)ziggurat_SHR3(anEngine);
96  iz=hz&127;
97  if((unsigned long)std::abs(hz)<kn[iz]) return (hz*wn[iz]);
98  }
99 }
100 
103 }
104 
105 double RandGaussZiggurat::operator()( double mean, double stdDev ) {
106  return ziggurat_RNOR(localEngine.get()) * stdDev + mean;
107 }
108 
109 void RandGaussZiggurat::shootArray( const int size, float* vect, float mean, float stdDev )
110 {
111  for (int i=0; i<size; ++i) {
112  vect[i] = shoot(mean,stdDev);
113  }
114 }
115 
116 void RandGaussZiggurat::shootArray( const int size, double* vect, double mean, double stdDev )
117 {
118  for (int i=0; i<size; ++i) {
119  vect[i] = shoot(mean,stdDev);
120  }
121 }
122 
123 void RandGaussZiggurat::shootArray( HepRandomEngine* anEngine, const int size, float* vect, float mean, float stdDev )
124 {
125  for (int i=0; i<size; ++i) {
126  vect[i] = shoot(anEngine,mean,stdDev);
127  }
128 }
129 
130 void RandGaussZiggurat::shootArray( HepRandomEngine* anEngine, const int size, double* vect, double mean, double stdDev )
131 {
132  for (int i=0; i<size; ++i) {
133  vect[i] = shoot(anEngine,mean,stdDev);
134  }
135 }
136 
137 void RandGaussZiggurat::fireArray( const int size, float* vect)
138 {
139  for (int i=0; i<size; ++i) {
140  vect[i] = fire( defaultMean, defaultStdDev );
141  }
142 }
143 
144 void RandGaussZiggurat::fireArray( const int size, double* vect)
145 {
146  for (int i=0; i<size; ++i) {
147  vect[i] = fire( defaultMean, defaultStdDev );
148  }
149 }
150 
151 void RandGaussZiggurat::fireArray( const int size, float* vect, float mean, float stdDev )
152 {
153  for (int i=0; i<size; ++i) {
154  vect[i] = fire( mean, stdDev );
155  }
156 }
157 
158 void RandGaussZiggurat::fireArray( const int size, double* vect, double mean, double stdDev )
159 {
160  for (int i=0; i<size; ++i) {
161  vect[i] = fire( mean, stdDev );
162  }
163 }
164 
165 std::ostream & RandGaussZiggurat::put ( std::ostream & os ) const {
166  int pr=os.precision(20);
167  os << " " << name() << "\n";
168  RandGauss::put(os);
169  os.precision(pr);
170  return os;
171 }
172 
173 std::istream & RandGaussZiggurat::get ( std::istream & is ) {
174  std::string inName;
175  is >> inName;
176  if (inName != name()) {
177  is.clear(std::ios::badbit | is.rdstate());
178  std::cerr << "Mismatch when expecting to read state of a "
179  << name() << " distribution\n"
180  << "Name found was " << inName
181  << "\nistream is left in the badbit state\n";
182  return is;
183  }
184  RandGauss::get(is);
185  return is;
186 }
187 
188 } // namespace CLHEP