Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/RandStudentT.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/RandStudentT.cc (Version 11.3.0) and /externals/clhep/src/RandStudentT.cc (Version 10.6.p1)


  1 // -*- C++ -*-                                      1 // -*- C++ -*-
  2 //                                                  2 //
  3 // -------------------------------------------      3 // -----------------------------------------------------------------------
  4 //                             HEP Random           4 //                             HEP Random
  5 //                          --- RandStudentT -      5 //                          --- RandStudentT ---
  6 //                      class implementation f      6 //                      class implementation file
  7 // -------------------------------------------      7 // -----------------------------------------------------------------------
  8                                                     8 
  9 // ===========================================      9 // =======================================================================
 10 // John Marraffino - Created: 12th May 1998        10 // John Marraffino - Created: 12th May 1998
 11 // G.Cosmo         - Fixed minor bug on inline     11 // G.Cosmo         - Fixed minor bug on inline definition for shoot()
 12 //                   methods : 20th Aug 1998       12 //                   methods : 20th Aug 1998
 13 // M Fischler      - put and get to/from strea     13 // M Fischler      - put and get to/from streams 12/13/04
 14 // M Fischler       - put/get to/from streams      14 // M Fischler       - put/get to/from streams uses pairs of ulongs when
 15 //      + storing doubles avoid problems with      15 //      + storing doubles avoid problems with precision 
 16 //      4/14/05                                    16 //      4/14/05
 17 // ===========================================     17 // =======================================================================
 18                                                    18 
 19 #include <float.h>                                 19 #include <float.h>
                                                   >>  20 #include <cmath>  // for std::log() std::exp()
 20 #include "CLHEP/Random/RandStudentT.h"             21 #include "CLHEP/Random/RandStudentT.h"
 21 #include "CLHEP/Random/DoubConv.h"                 22 #include "CLHEP/Random/DoubConv.h"
 22                                                    23 
 23 #include <cmath>  // for std::log() std::exp() << 
 24 #include <iostream>                            << 
 25 #include <string>                              << 
 26 #include <vector>                              << 
 27                                                << 
 28 namespace CLHEP {                                  24 namespace CLHEP {
 29                                                    25 
 30 std::string RandStudentT::name() const {return     26 std::string RandStudentT::name() const {return "RandStudentT";}
 31 HepRandomEngine & RandStudentT::engine() {retu     27 HepRandomEngine & RandStudentT::engine() {return *localEngine;}
 32                                                    28 
 33 RandStudentT::~RandStudentT()                      29 RandStudentT::~RandStudentT()
 34 {                                                  30 {
 35 }                                                  31 }
 36                                                    32 
 37 double RandStudentT::operator()() {                33 double RandStudentT::operator()() {
 38   return fire( defaultA );                         34   return fire( defaultA );
 39 }                                                  35 }
 40                                                    36 
 41 double RandStudentT::operator()( double a ) {      37 double RandStudentT::operator()( double a ) {
 42   return fire( a );                                38   return fire( a );
 43 }                                                  39 }
 44                                                    40 
 45 double RandStudentT::shoot( double a ) {           41 double RandStudentT::shoot( double a ) {
 46 /*********************************************     42 /******************************************************************
 47  *                                                 43  *                                                                *
 48  *           Student-t Distribution - Polar Me     44  *           Student-t Distribution - Polar Method                *
 49  *                                                 45  *                                                                *
 50  *********************************************     46  ******************************************************************
 51  *                                                 47  *                                                                *
 52  * The polar method of Box/Muller for generati     48  * The polar method of Box/Muller for generating Normal variates  *
 53  * is adapted to the Student-t distribution. T     49  * is adapted to the Student-t distribution. The two generated    *
 54  * variates are not independent and the expect     50  * variates are not independent and the expected no. of uniforms  *
 55  * per variate is 2.5464.                          51  * per variate is 2.5464.                                         *
 56  *                                                 52  *                                                                *
 57  *********************************************     53  ******************************************************************
 58  *                                                 54  *                                                                *
 59  * FUNCTION :   - tpol  samples a random numbe     55  * FUNCTION :   - tpol  samples a random number from the          *
 60  *                Student-t distribution with      56  *                Student-t distribution with parameter a > 0.    *
 61  * REFERENCE :  - R.W. Bailey (1994): Polar ge     57  * REFERENCE :  - R.W. Bailey (1994): Polar generation of random  *
 62  *                variates with the t-distribu     58  *                variates with the t-distribution, Mathematics   *
 63  *                of Computation 62, 779-781.      59  *                of Computation 62, 779-781.                     *
 64  * SUBPROGRAM : -  ... (0,1)-Uniform generator     60  * SUBPROGRAM : -  ... (0,1)-Uniform generator                    *
 65  *                                                 61  *                                                                *
 66  *                                                 62  *                                                                *
 67  * Implemented by F. Niederl, 1994                 63  * Implemented by F. Niederl, 1994                                *
 68  *********************************************     64  ******************************************************************/
 69                                                    65 
 70  double u,v,w;                                     66  double u,v,w;
 71                                                    67 
 72 // Check for valid input value                     68 // Check for valid input value
 73                                                    69 
 74  if( a < 0.0) return (DBL_MAX);                    70  if( a < 0.0) return (DBL_MAX);
 75                                                    71 
 76  do                                                72  do
 77  {                                                 73  {
 78    u = 2.0 * HepRandom::getTheEngine()->flat()     74    u = 2.0 * HepRandom::getTheEngine()->flat() - 1.0;
 79    v = 2.0 * HepRandom::getTheEngine()->flat()     75    v = 2.0 * HepRandom::getTheEngine()->flat() - 1.0;
 80  }                                                 76  }
 81  while ((w = u * u + v * v) > 1.0);                77  while ((w = u * u + v * v) > 1.0);
 82                                                    78 
 83  return(u * std::sqrt( a * ( std::exp(- 2.0 /      79  return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
 84 }                                                  80 }
 85                                                    81 
 86 void RandStudentT::shootArray( const int size,     82 void RandStudentT::shootArray( const int size, double* vect,
 87                             double a )             83                             double a )
 88 {                                                  84 {
 89   for( double* v = vect; v != vect + size; ++v     85   for( double* v = vect; v != vect + size; ++v )
 90     *v = shoot(a);                                 86     *v = shoot(a);
 91 }                                                  87 }
 92                                                    88 
 93 void RandStudentT::shootArray( HepRandomEngine     89 void RandStudentT::shootArray( HepRandomEngine* anEngine,
 94                             const int size, do     90                             const int size, double* vect,
 95                             double a )             91                             double a )
 96 {                                                  92 {
 97   for( double* v = vect; v != vect + size; ++v     93   for( double* v = vect; v != vect + size; ++v )
 98     *v = shoot(anEngine,a);                        94     *v = shoot(anEngine,a);
 99 }                                                  95 }
100                                                    96 
101 double RandStudentT::fire( double a ) {            97 double RandStudentT::fire( double a ) {
102                                                    98 
103  double u,v,w;                                     99  double u,v,w;
104                                                   100 
105  do                                               101  do
106  {                                                102  {
107    u = 2.0 * localEngine->flat() - 1.0;           103    u = 2.0 * localEngine->flat() - 1.0;
108    v = 2.0 * localEngine->flat() - 1.0;           104    v = 2.0 * localEngine->flat() - 1.0;
109  }                                                105  }
110  while ((w = u * u + v * v) > 1.0);               106  while ((w = u * u + v * v) > 1.0);
111                                                   107 
112  return(u * std::sqrt( a * ( std::exp(- 2.0 /     108  return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
113 }                                                 109 }
114                                                   110 
115 void RandStudentT::fireArray( const int size,     111 void RandStudentT::fireArray( const int size, double* vect)
116 {                                                 112 {
117   for( double* v = vect; v != vect + size; ++v    113   for( double* v = vect; v != vect + size; ++v )
118     *v = fire(defaultA);                          114     *v = fire(defaultA);
119 }                                                 115 }
120                                                   116 
121 void RandStudentT::fireArray( const int size,     117 void RandStudentT::fireArray( const int size, double* vect,
122                            double a )             118                            double a )
123 {                                                 119 {
124   for( double* v = vect; v != vect + size; ++v    120   for( double* v = vect; v != vect + size; ++v )
125     *v = fire(a);                                 121     *v = fire(a);
126 }                                                 122 }
127                                                   123 
128 double RandStudentT::shoot( HepRandomEngine *a    124 double RandStudentT::shoot( HepRandomEngine *anEngine, double a ) {
129                                                   125 
130  double u,v,w;                                    126  double u,v,w;
131                                                   127 
132  do                                               128  do
133  {                                                129  {
134    u = 2.0 * anEngine->flat() - 1.0;              130    u = 2.0 * anEngine->flat() - 1.0;
135    v = 2.0 * anEngine->flat() - 1.0;              131    v = 2.0 * anEngine->flat() - 1.0;
136  }                                                132  }
137  while ((w = u * u + v * v) > 1.0);               133  while ((w = u * u + v * v) > 1.0);
138                                                   134 
139  return(u * std::sqrt( a * ( std::exp(- 2.0 /     135  return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
140 }                                                 136 }
141                                                   137 
142 std::ostream & RandStudentT::put ( std::ostrea    138 std::ostream & RandStudentT::put ( std::ostream & os ) const {
143   long pr=os.precision(20);                    << 139   int pr=os.precision(20);
144   std::vector<unsigned long> t(2);                140   std::vector<unsigned long> t(2);
145   os << " " << name() << "\n";                    141   os << " " << name() << "\n";
146   os << "Uvec" << "\n";                           142   os << "Uvec" << "\n";
147   t = DoubConv::dto2longs(defaultA);              143   t = DoubConv::dto2longs(defaultA);
148   os << defaultA << " " << t[0] << " " << t[1]    144   os << defaultA << " " << t[0] << " " << t[1] << "\n";
149   os.precision(pr);                               145   os.precision(pr);
150   return os;                                      146   return os;
151 }                                                 147 }
152                                                   148 
153 std::istream & RandStudentT::get ( std::istrea    149 std::istream & RandStudentT::get ( std::istream & is ) {
154   std::string inName;                             150   std::string inName;
155   is >> inName;                                   151   is >> inName;
156   if (inName != name()) {                         152   if (inName != name()) {
157     is.clear(std::ios::badbit | is.rdstate());    153     is.clear(std::ios::badbit | is.rdstate());
158     std::cerr << "Mismatch when expecting to r    154     std::cerr << "Mismatch when expecting to read state of a "
159             << name() << " distribution\n"        155             << name() << " distribution\n"
160         << "Name found was " << inName            156         << "Name found was " << inName
161         << "\nistream is left in the badbit st    157         << "\nistream is left in the badbit state\n";
162     return is;                                    158     return is;
163   }                                               159   }
164   if (possibleKeywordInput(is, "Uvec", default    160   if (possibleKeywordInput(is, "Uvec", defaultA)) {
165     std::vector<unsigned long> t(2);              161     std::vector<unsigned long> t(2);
166     is >> defaultA >> t[0] >> t[1]; defaultA =    162     is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t); 
167     return is;                                    163     return is;
168   }                                               164   }
169   // is >> defaultA encompassed by possibleKey    165   // is >> defaultA encompassed by possibleKeywordInput
170   return is;                                      166   return is;
171 }                                                 167 }
172                                                   168 
173 }  // namespace CLHEP                             169 }  // namespace CLHEP
174                                                   170