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  }
191  return is;
192 }
193 
194 } // namespace CLHEP
CLHEP::RandGaussQ::engine
HepRandomEngine & engine()
Definition: RandGaussQ.cc:24
CLHEP::RandGauss::engine
HepRandomEngine & engine()
Definition: RandGauss.cc:44
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
CLHEP::RandGaussQ::fire
double fire()
Table1size
#define Table1size
Definition: RandGaussQ.cc:77
is
HepRotation and so forth isNear() norm2() rectify() static Rotation row1 row4(To avoid bloat in the code pulled in for programs which don 't use all these features, we split the implementation .cc files. Only isNear() goes into the original Rotation.cc) --------------------------------------- HepAxisAngle and HepEulerAngles classes --------------------------------------- These classes are very useful and simple structures for holding the result of a nice intuituve decomposition of a rotation there is no longer much content in the distinct ZOOM PhysicsVectors library The only content left in the library is the object files representing the various Exception objects When we build the CLHEP classes for the ZOOM we will set up so as to use ZOOM SpaceVector is(but we can disable namespace usage and most of our users do so at this point). What I do is leave Hep3Vector in the global namespace
CLHEP::RandGaussQ::shoot
static double shoot()
Table0offset
#define Table0offset
Definition: RandGaussQ.cc:85
CLHEP::RandGauss::get
std::istream & get(std::istream &is)
Definition: RandGauss.cc:262
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
Table0size
#define Table0size
Definition: RandGaussQ.cc:76
CLHEP::RandGauss::defaultStdDev
double defaultStdDev
Definition: Matrix/CLHEP/Random/RandGauss.h:153
CLHEP
Definition: ClhepVersion.h:13
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
CLHEP::RandGauss::defaultMean
double defaultMean
Definition: Matrix/CLHEP/Random/RandGauss.h:152
Table0step
#define Table0step
Definition: RandGaussQ.cc:80
Table0scale
#define Table0scale
Definition: RandGaussQ.cc:83
CLHEP::RandGaussQ::shootArray
static void shootArray(const int size, double *vect, double mean=0.0, double stdDev=1.0)
Definition: RandGaussQ.cc:37
CLHEP::RandGaussQ::transformQuick
static double transformQuick(double r)
Definition: RandGaussQ.cc:95
Table1offset
#define Table1offset
Definition: RandGaussQ.cc:86
CLHEP::RandGaussQ::transformSmall
static double transformSmall(double r)
Definition: RandGaussQ.cc:130
CLHEP::RandGauss::put
std::ostream & put(std::ostream &os) const
Definition: RandGauss.cc:237
CLHEP::RandGaussQ::fireArray
void fireArray(const int size, double *vect)
Definition: RandGaussQ.cc:52
i
long i
Definition: JamesRandomSeeding.txt:27
TableSize
#define TableSize
Definition: RandGaussQ.cc:78
CLHEP::RandGaussQ::name
std::string name() const
Definition: RandGaussQ.cc:23
Table1step
#define Table1step
Definition: RandGaussQ.cc:81
CLHEP::RandGauss::localEngine
shared_ptr< HepRandomEngine > localEngine
Definition: Matrix/CLHEP/Random/RandGauss.h:155
CLHEP::RandGaussQ::~RandGaussQ
virtual ~RandGaussQ()
Definition: RandGaussQ.cc:26
CLHEP::RandGaussQ::get
std::istream & get(std::istream &is)
Definition: RandGaussQ.cc:179
CLHEP::RandGaussQ::operator()
virtual double operator()()
Definition: RandGaussQ.cc:29
CLHEP::RandGaussQ::put
std::ostream & put(std::ostream &os) const
Definition: RandGaussQ.cc:171