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