Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/RandExpZiggurat.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 #include "CLHEP/Random/DoubConv.h"
  2 
  3 #include "CLHEP/Random/RandExpZiggurat.h"
  4 
  5 #include <cmath>  // for std::log()
  6 #include <iostream>
  7 #include <vector>
  8 
  9 namespace CLHEP {
 10 
 11 CLHEP_THREAD_LOCAL unsigned long RandExpZiggurat::kn[128], RandExpZiggurat::ke[256];
 12 CLHEP_THREAD_LOCAL float RandExpZiggurat::wn[128],RandExpZiggurat::fn[128],RandExpZiggurat::we[256],RandExpZiggurat::fe[256];
 13 CLHEP_THREAD_LOCAL bool RandExpZiggurat::ziggurat_is_init = false;
 14 
 15 std::string RandExpZiggurat::name() const {return "RandExpZiggurat";}
 16 
 17 HepRandomEngine & RandExpZiggurat::engine() {return *localEngine;}
 18 
 19 RandExpZiggurat::~RandExpZiggurat() {
 20 }
 21 
 22 RandExpZiggurat::RandExpZiggurat(const RandExpZiggurat& right) : HepRandom(right),defaultMean(right.defaultMean)
 23 {
 24 }
 25 
 26 double RandExpZiggurat::operator()()
 27 {
 28   return fire( defaultMean );
 29 }
 30 
 31 void RandExpZiggurat::shootArray( const int size, float* vect, float mean )
 32 {
 33    for (int i=0; i<size; ++i) vect[i] = shoot(mean);
 34 }
 35 
 36 void RandExpZiggurat::shootArray( const int size, double* vect, double mean )
 37 {
 38    for (int i=0; i<size; ++i) vect[i] = shoot(mean);
 39 }
 40 
 41 void RandExpZiggurat::shootArray(HepRandomEngine* anEngine, const int size, float* vect, float mean )
 42 {
 43    for (int i=0; i<size; ++i) vect[i] = shoot(anEngine, mean);
 44 }
 45 
 46 void RandExpZiggurat::shootArray(HepRandomEngine* anEngine, const int size, double* vect, double mean )
 47 {
 48    for (int i=0; i<size; ++i) vect[i] = shoot(anEngine, mean);
 49 }
 50 
 51 void RandExpZiggurat::fireArray( const int size, float* vect)
 52 {
 53    for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
 54 }
 55 
 56 void RandExpZiggurat::fireArray( const int size, double* vect)
 57 {
 58    for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
 59 }
 60 
 61 void RandExpZiggurat::fireArray( const int size, float* vect, float mean )
 62 {
 63    for (int i=0; i<size; ++i) vect[i] = fire( mean );
 64 }
 65 
 66 void RandExpZiggurat::fireArray( const int size, double* vect, double mean )
 67 {
 68    for (int i=0; i<size; ++i) vect[i] = fire( mean );
 69 }
 70 
 71 std::ostream & RandExpZiggurat::put ( std::ostream & os ) const {
 72   int pr=os.precision(20);
 73   std::vector<unsigned long> t(2);
 74   os << " " << name() << "\n";
 75   os << "Uvec" << "\n";
 76   t = DoubConv::dto2longs(defaultMean);
 77   os << defaultMean << " " << t[0] << " " << t[1] << "\n";
 78   os.precision(pr);
 79   return os;
 80 #ifdef REMOVED
 81   int pr=os.precision(20);
 82   os << " " << name() << "\n";
 83   os << defaultMean << "\n";
 84   os.precision(pr);
 85   return os;
 86 #endif
 87 }
 88 
 89 std::istream & RandExpZiggurat::get ( std::istream & is ) {
 90   std::string inName;
 91   is >> inName;
 92   if (inName != name()) {
 93     is.clear(std::ios::badbit | is.rdstate());
 94     std::cerr << "Mismatch when expecting to read state of a "
 95             << name() << " distribution\n"
 96         << "Name found was " << inName
 97         << "\nistream is left in the badbit state\n";
 98     return is;
 99   }
100   if (possibleKeywordInput(is, "Uvec", defaultMean)) {
101     std::vector<unsigned long> t(2);
102     is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t); 
103     return is;
104   }
105   // is >> defaultMean encompassed by possibleKeywordInput
106   return is;
107 }
108 
109 
110 float RandExpZiggurat::ziggurat_efix(unsigned long jz,HepRandomEngine* anEngine)
111 { 
112   if(!ziggurat_is_init) ziggurat_init();
113 
114   unsigned long iz=jz&255;
115   
116   float x;
117   for(;;)
118   {  
119     if(iz==0) return (7.69711-std::log(ziggurat_UNI(anEngine)));          /* iz==0 */
120     x=jz*we[iz]; 
121     if( fe[iz]+ziggurat_UNI(anEngine)*(fe[iz-1]-fe[iz]) < std::exp(-x) ) return (x);
122 
123     /* initiate, try to exit for(;;) loop */
124     jz=ziggurat_SHR3(anEngine);
125     iz=(jz&255);
126     if(jz<ke[iz]) return (jz*we[iz]);
127   }
128 }
129 
130 bool RandExpZiggurat::ziggurat_init()
131 {  
132   const double rzm1 = 2147483648.0, rzm2 = 4294967296.;
133   double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
134   double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
135   int i;
136 
137 /* Set up tables for RNOR */
138   q=vn/std::exp(-.5*dn*dn);
139   kn[0]=(unsigned long)((dn/q)*rzm1);
140   kn[1]=0;
141 
142   wn[0]=q/rzm1;
143   wn[127]=dn/rzm1;
144 
145   fn[0]=1.;
146   fn[127]=std::exp(-.5*dn*dn);
147 
148   for(i=126;i>=1;i--) {
149     dn=std::sqrt(-2.*std::log(vn/dn+std::exp(-.5*dn*dn)));
150     kn[i+1]=(unsigned long)((dn/tn)*rzm1);
151     tn=dn;
152     fn[i]=std::exp(-.5*dn*dn);
153     wn[i]=dn/rzm1;
154   }
155 
156 /* Set up tables for REXP */
157   q = ve/std::exp(-de);
158   ke[0]=(unsigned long)((de/q)*rzm2);
159   ke[1]=0;
160 
161   we[0]=q/rzm2;
162   we[255]=de/rzm2;
163 
164   fe[0]=1.;
165   fe[255]=std::exp(-de);
166 
167   for(i=254;i>=1;i--) {
168     de=-std::log(ve/de+std::exp(-de));
169     ke[i+1]= (unsigned long)((de/te)*rzm2);
170     te=de;
171     fe[i]=std::exp(-de);
172     we[i]=de/rzm2;
173   }
174   ziggurat_is_init=true;
175   return true;
176 }
177 
178 }  // namespace CLHEP
179