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