ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RandGaussQ.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RandGaussQ.cc
1 // -*- C++ -*-
2 //
3 // -----------------------------------------------------------------------
4 // HEP Random
5 // --- RandGaussQ ---
6 // class implementation file
7 // -----------------------------------------------------------------------
8 
9 // =======================================================================
10 // M Fischler - Created 24 Jan 2000
11 // M Fischler - put and get to/from streams 12/13/04
12 // =======================================================================
13 
16 #include <iostream>
17 #include <cmath> // for std::log()
18 
19 namespace CLHEP {
20 
21 std::string RandGaussQ::name() const {return "RandGaussQ";}
23 
25 }
26 
29 }
30 
31 double RandGaussQ::operator()( double mean, double stdDev ) {
32  return transformQuick(localEngine->flat()) * stdDev + mean;
33 }
34 
35 void RandGaussQ::shootArray( const int size, double* vect,
36  double mean, double stdDev )
37 {
38  for( double* v = vect; v != vect + size; ++v )
39  *v = shoot(mean,stdDev);
40 }
41 
43  const int size, double* vect,
44  double mean, double stdDev )
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, double* vect)
51 {
52  for( double* v = vect; v != vect + size; ++v )
54 }
55 
56 void RandGaussQ::fireArray( const int size, double* vect,
57  double mean, double stdDev )
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) and quickTransform(r)
65 //
66 
67 // Since all these are this is static to this compilation unit only, the
68 // info is establised a priori and not at each invocation.
69 
70 // The main data is of course the gaussQTables table; the rest is all
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 in a file ---
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 negative number of
93  // sigmas. For r > 0 we will multiply by
94  // sign = -1 to return a positive number.
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 Table1size
105  if (index == Table1size) return 0.0;
106  dx = (Table1size<<1) * r - index; // fraction of way to next bin
107  index += Table1offset-1;
108  } else if ( r > Table0step ) {
109  double rr = r * Table0scale;
110  index = int(Table0size * rr); // 1 to Table0size
111  dx = Table0size * rr - index; // fraction of way to next bin
112  index += Table0offset-1;
113  } else { // r <= Table0step - not in tables
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-dx) ));
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 1*3 1*3*5
129  // ------------ * (1 - ---- + ---- - ----- + ... )
130  // v*sqrt(2*pi) v**2 v**4 v**6
131 
132  // 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 value of r is meaningful only
134  // to an absolute error of 1E-16 (double precision accuracy for a number
135  // which on the high side could be of the form 1-epsilon), computing
136  // v to more than 3-4 digits of accuracy is suspect; however, to ensure
137  // smoothness with the table generator (which uses quite a few terms) we
138  // 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 
141  // This routine is called less than one time in a million firings, so
142  // speed is of no concern. As a matter of technique, we terminate the
143  // iterations in case they would be infinite, but this should not happen.
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*vn2*vn2*vn2*vn2;
152  s1 += 11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2;
153  s1 += -9*7*5*3 * vn2*vn2*vn2*vn2*vn2;
154  s1 += 7*5*3 * vn2*vn2*vn2*vn2;
155  s1 += -5*3 * vn2*vn2*vn2;
156  s1 += 3 * vn2*vn2 - vn2 + 1.0;
157  v = std::sqrt ( 2.0 * std::log ( s1 / (r*guess*std::sqrt(CLHEP::twopi)) ) );
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 & 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 & RandGaussQ::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