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