CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

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