Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/RandGaussQ.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /externals/clhep/src/RandGaussQ.cc (Version 11.3.0) and /externals/clhep/src/RandGaussQ.cc (Version 9.6.p3)


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