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