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