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