| Geant4 Cross Reference |
>> 1 // $Id:$
1 // -*- C++ -*- 2 // -*- C++ -*-
2 // 3 //
3 // ------------------------------------------- 4 // -----------------------------------------------------------------------
4 // HEP Random 5 // HEP Random
5 // --- RandGaussQ --- 6 // --- RandGaussQ ---
6 // class implementation f 7 // class implementation file
7 // ------------------------------------------- 8 // -----------------------------------------------------------------------
8 9
9 // =========================================== 10 // =======================================================================
10 // M Fischler - Created 24 Jan 2000 11 // M Fischler - Created 24 Jan 2000
11 // M Fischler - put and get to/from stream 12 // M Fischler - put and get to/from streams 12/13/04
12 // =========================================== 13 // =======================================================================
13 14
14 #include "CLHEP/Random/RandGaussQ.h" 15 #include "CLHEP/Random/RandGaussQ.h"
15 #include "CLHEP/Units/PhysicalConstants.h" 16 #include "CLHEP/Units/PhysicalConstants.h"
16 #include <iostream> 17 #include <iostream>
17 #include <cmath> // for std::log() 18 #include <cmath> // for std::log()
18 19
19 namespace CLHEP { 20 namespace CLHEP {
20 21
21 std::string RandGaussQ::name() const {return " 22 std::string RandGaussQ::name() const {return "RandGaussQ";}
22 HepRandomEngine & RandGaussQ::engine() {return 23 HepRandomEngine & RandGaussQ::engine() {return RandGauss::engine();}
23 24
24 RandGaussQ::~RandGaussQ() { 25 RandGaussQ::~RandGaussQ() {
25 } 26 }
26 27
27 double RandGaussQ::operator()() { 28 double RandGaussQ::operator()() {
28 return transformQuick(localEngine->flat()) * 29 return transformQuick(localEngine->flat()) * defaultStdDev + defaultMean;
29 } 30 }
30 31
31 double RandGaussQ::operator()( double mean, do 32 double RandGaussQ::operator()( double mean, double stdDev ) {
32 return transformQuick(localEngine->flat()) * 33 return transformQuick(localEngine->flat()) * stdDev + mean;
33 } 34 }
34 35
35 void RandGaussQ::shootArray( const int size, d 36 void RandGaussQ::shootArray( const int size, double* vect,
36 double mean, doubl 37 double mean, double stdDev )
37 { 38 {
38 for( double* v = vect; v != vect + size; ++v 39 for( double* v = vect; v != vect + size; ++v )
39 *v = shoot(mean,stdDev); 40 *v = shoot(mean,stdDev);
40 } 41 }
41 42
42 void RandGaussQ::shootArray( HepRandomEngine* 43 void RandGaussQ::shootArray( HepRandomEngine* anEngine,
43 const int size, do 44 const int size, double* vect,
44 double mean, doubl 45 double mean, double stdDev )
45 { 46 {
46 for( double* v = vect; v != vect + size; ++v 47 for( double* v = vect; v != vect + size; ++v )
47 *v = shoot(anEngine,mean,stdDev); 48 *v = shoot(anEngine,mean,stdDev);
48 } 49 }
49 50
50 void RandGaussQ::fireArray( const int size, do 51 void RandGaussQ::fireArray( const int size, double* vect)
51 { 52 {
52 for( double* v = vect; v != vect + size; ++v 53 for( double* v = vect; v != vect + size; ++v )
53 *v = fire( defaultMean, defaultStdDev ); 54 *v = fire( defaultMean, defaultStdDev );
54 } 55 }
55 56
56 void RandGaussQ::fireArray( const int size, do 57 void RandGaussQ::fireArray( const int size, double* vect,
57 double mean, double 58 double mean, double stdDev )
58 { 59 {
59 for( double* v = vect; v != vect + size; ++v 60 for( double* v = vect; v != vect + size; ++v )
60 *v = fire( mean, stdDev ); 61 *v = fire( mean, stdDev );
61 } 62 }
62 63
63 // 64 //
64 // Table of errInts, for use with transform(r) 65 // Table of errInts, for use with transform(r) and quickTransform(r)
65 // 66 //
66 67
67 // Since all these are this is static to this 68 // Since all these are this is static to this compilation unit only, the
68 // info is establised a priori and not at each 69 // info is establised a priori and not at each invocation.
69 70
70 // The main data is of course the gaussQTables 71 // The main data is of course the gaussQTables table; the rest is all
71 // bookkeeping to know what the tables mean. 72 // bookkeeping to know what the tables mean.
72 73
73 #define Table0size 250 74 #define Table0size 250
74 #define Table1size 1000 75 #define Table1size 1000
75 #define TableSize (Table0size+Table1size) 76 #define TableSize (Table0size+Table1size)
76 77
77 #define Table0step (2.0E-6) 78 #define Table0step (2.0E-6)
78 #define Table1step (5.0E-4) 79 #define Table1step (5.0E-4)
79 80
80 #define Table0scale (1.0/Table1step) 81 #define Table0scale (1.0/Table1step)
81 82
82 #define Table0offset 0 83 #define Table0offset 0
83 #define Table1offset (Table0size) 84 #define Table1offset (Table0size)
84 85
85 // Here comes the big (5K bytes) table, kept 86 // Here comes the big (5K bytes) table, kept in a file ---
86 87
87 static const float gaussTables [TableSize] = { 88 static const float gaussTables [TableSize] = {
88 #include "CLHEP/Random/gaussQTables.cdat" 89 #include "CLHEP/Random/gaussQTables.cdat"
89 }; 90 };
90 91
91 double RandGaussQ::transformQuick (double r) { 92 double RandGaussQ::transformQuick (double r) {
92 double sign = +1.0; // We always compute a n 93 double sign = +1.0; // We always compute a negative number of
93 // sigmas. For r > 0 we will multiply 94 // sigmas. For r > 0 we will multiply by
94 // sign = -1 to return a positive numb 95 // sign = -1 to return a positive number.
95 if ( r > .5 ) { 96 if ( r > .5 ) {
96 r = 1-r; 97 r = 1-r;
97 sign = -1.0; 98 sign = -1.0;
98 } 99 }
99 100
100 int index; 101 int index;
101 double dx; 102 double dx;
102 103
103 if ( r >= Table1step ) { 104 if ( r >= Table1step ) {
104 index = int((Table1size<<1) * r); // 1 to 105 index = int((Table1size<<1) * r); // 1 to Table1size
105 if (index == Table1size) return 0.0; 106 if (index == Table1size) return 0.0;
106 dx = (Table1size<<1) * r - index; // f 107 dx = (Table1size<<1) * r - index; // fraction of way to next bin
107 index += Table1offset-1; 108 index += Table1offset-1;
108 } else if ( r > Table0step ) { 109 } else if ( r > Table0step ) {
109 double rr = r * Table0scale; 110 double rr = r * Table0scale;
110 index = int(Table0size * rr); // 1 to Ta 111 index = int(Table0size * rr); // 1 to Table0size
111 dx = Table0size * rr - index; // fract 112 dx = Table0size * rr - index; // fraction of way to next bin
112 index += Table0offset-1; 113 index += Table0offset-1;
113 } else { // r <= Table0step - not 114 } else { // r <= Table0step - not in tables
114 return sign*transformSmall(r); 115 return sign*transformSmall(r);
115 } 116 }
116 117
117 double y0 = gaussTables [index++]; 118 double y0 = gaussTables [index++];
118 double y1 = gaussTables [index]; 119 double y1 = gaussTables [index];
119 120
120 return (float) (sign * ( y1 * dx + y0 * (1.0 121 return (float) (sign * ( y1 * dx + y0 * (1.0-dx) ));
121 122
122 } // transformQuick() 123 } // transformQuick()
123 124
124 double RandGaussQ::transformSmall (double r) { 125 double RandGaussQ::transformSmall (double r) {
125 126
126 // Solve for -v in the asymtotic formula 127 // Solve for -v in the asymtotic formula
127 // 128 //
128 // errInt (-v) = exp(-v*v/2) 1 129 // errInt (-v) = exp(-v*v/2) 1 1*3 1*3*5
129 // ------------ * (1 - ---- + ---- - - 130 // ------------ * (1 - ---- + ---- - ----- + ... )
130 // v*sqrt(2*pi) v**2 v**4 v 131 // v*sqrt(2*pi) v**2 v**4 v**6
131 132
132 // The value of r (=errInt(-v)) supplied is 133 // The value of r (=errInt(-v)) supplied is going to less than 2.0E-13,
133 // which is such that v < -7.25. Since the 134 // which is such that v < -7.25. Since the value of r is meaningful only
134 // to an absolute error of 1E-16 (double pre 135 // to an absolute error of 1E-16 (double precision accuracy for a number
135 // which on the high side could be of the fo 136 // which on the high side could be of the form 1-epsilon), computing
136 // v to more than 3-4 digits of accuracy is 137 // v to more than 3-4 digits of accuracy is suspect; however, to ensure
137 // smoothness with the table generator (whic 138 // smoothness with the table generator (which uses quite a few terms) we
138 // also use terms up to 1*3*5* ... *13/v**14 139 // also use terms up to 1*3*5* ... *13/v**14, and insist on accuracy of
139 // solution at the level of 1.0e-7. 140 // solution at the level of 1.0e-7.
140 141
141 // This routine is called less than one time 142 // This routine is called less than one time in a million firings, so
142 // speed is of no concern. As a matter of t 143 // speed is of no concern. As a matter of technique, we terminate the
143 // iterations in case they would be infinite 144 // iterations in case they would be infinite, but this should not happen.
144 145
145 double eps = 1.0e-7; 146 double eps = 1.0e-7;
146 double guess = 7.5; 147 double guess = 7.5;
147 double v; 148 double v;
148 149
149 for ( int i = 1; i < 50; i++ ) { 150 for ( int i = 1; i < 50; i++ ) {
150 double vn2 = 1.0/(guess*guess); 151 double vn2 = 1.0/(guess*guess);
151 double s1 = -13*11*9*7*5*3 * vn2*vn2*vn2*v 152 double s1 = -13*11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2*vn2;
152 s1 += 11*9*7*5*3 * vn2*vn2*vn2* 153 s1 += 11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2;
153 s1 += -9*7*5*3 * vn2*vn2*vn2* 154 s1 += -9*7*5*3 * vn2*vn2*vn2*vn2*vn2;
154 s1 += 7*5*3 * vn2*vn2*vn2* 155 s1 += 7*5*3 * vn2*vn2*vn2*vn2;
155 s1 += -5*3 * vn2*vn2*vn2; 156 s1 += -5*3 * vn2*vn2*vn2;
156 s1 += 3 * vn2*vn2 - 157 s1 += 3 * vn2*vn2 - vn2 + 1.0;
157 v = std::sqrt ( 2.0 * std::log ( s1 / (r*g 158 v = std::sqrt ( 2.0 * std::log ( s1 / (r*guess*std::sqrt(CLHEP::twopi)) ) );
158 if ( std::fabs(v-guess) < eps ) break; 159 if ( std::fabs(v-guess) < eps ) break;
159 guess = v; 160 guess = v;
160 } 161 }
161 return -v; 162 return -v;
162 163
163 } // transformSmall() 164 } // transformSmall()
164 165
165 std::ostream & RandGaussQ::put ( std::ostream 166 std::ostream & RandGaussQ::put ( std::ostream & os ) const {
166 long pr=os.precision(20); << 167 int pr=os.precision(20);
167 os << " " << name() << "\n"; 168 os << " " << name() << "\n";
168 RandGauss::put(os); 169 RandGauss::put(os);
169 os.precision(pr); 170 os.precision(pr);
170 return os; 171 return os;
171 } 172 }
172 173
173 std::istream & RandGaussQ::get ( std::istream 174 std::istream & RandGaussQ::get ( std::istream & is ) {
174 std::string inName; 175 std::string inName;
175 is >> inName; 176 is >> inName;
176 if (inName != name()) { 177 if (inName != name()) {
177 is.clear(std::ios::badbit | is.rdstate()); 178 is.clear(std::ios::badbit | is.rdstate());
178 std::cerr << "Mismatch when expecting to r 179 std::cerr << "Mismatch when expecting to read state of a "
179 << name() << " distribution\n" 180 << name() << " distribution\n"
180 << "Name found was " << inName 181 << "Name found was " << inName
181 << "\nistream is left in the badbit st 182 << "\nistream is left in the badbit state\n";
182 return is; 183 return is;
183 } 184 }
184 RandGauss::get(is); 185 RandGauss::get(is);
185 return is; 186 return is;
186 } 187 }
187 188
188 } // namespace CLHEP 189 } // namespace CLHEP
189 190