| Geant4 Cross Reference |
>> 1 // $Id:$
1 // -*- C++ -*- 2 // -*- C++ -*-
2 // 3 //
3 // ------------------------------------------- 4 // -----------------------------------------------------------------------
4 // HEP Random 5 // HEP Random
5 // --- RandChiSquare 6 // --- RandChiSquare ---
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 // M Fischler - put and get to/from stream 12 // M Fischler - put and get to/from streams 12/10/04
12 // M Fischler - put/get to/from streams 13 // M Fischler - put/get to/from streams uses pairs of ulongs when
13 // + storing doubles avoid problems with 14 // + storing doubles avoid problems with precision
14 // 4/14/05 15 // 4/14/05
15 // =========================================== 16 // =======================================================================
16 17
17 #include "CLHEP/Random/RandChiSquare.h" 18 #include "CLHEP/Random/RandChiSquare.h"
18 #include "CLHEP/Random/DoubConv.h" 19 #include "CLHEP/Random/DoubConv.h"
19 #include "CLHEP/Utility/thread_local.h" << 20 #include <cmath> // for log()
20 #include <cmath> // for std::log() <<
21 #include <iostream> <<
22 #include <string> <<
23 #include <vector> <<
24 21
25 namespace CLHEP { 22 namespace CLHEP {
26 23
27 std::string RandChiSquare::name() const {retur 24 std::string RandChiSquare::name() const {return "RandChiSquare";}
28 HepRandomEngine & RandChiSquare::engine() {ret 25 HepRandomEngine & RandChiSquare::engine() {return *localEngine;}
29 26
30 RandChiSquare::~RandChiSquare() { 27 RandChiSquare::~RandChiSquare() {
31 } 28 }
32 29
33 double RandChiSquare::shoot( HepRandomEngine * 30 double RandChiSquare::shoot( HepRandomEngine *anEngine, double a ) {
34 return genChiSquare( anEngine, a ); 31 return genChiSquare( anEngine, a );
35 } 32 }
36 33
37 double RandChiSquare::shoot( double a ) { 34 double RandChiSquare::shoot( double a ) {
38 HepRandomEngine *anEngine = HepRandom::getTh 35 HepRandomEngine *anEngine = HepRandom::getTheEngine();
39 return genChiSquare( anEngine, a ); 36 return genChiSquare( anEngine, a );
40 } 37 }
41 38
42 double RandChiSquare::fire( double a ) { 39 double RandChiSquare::fire( double a ) {
43 return genChiSquare( localEngine.get(), a ); 40 return genChiSquare( localEngine.get(), a );
44 } 41 }
45 42
46 void RandChiSquare::shootArray( const int size 43 void RandChiSquare::shootArray( const int size, double* vect,
47 double a ) { 44 double a ) {
48 for( double* v = vect; v != vect+size; ++v ) 45 for( double* v = vect; v != vect+size; ++v )
49 *v = shoot(a); 46 *v = shoot(a);
50 } 47 }
51 48
52 void RandChiSquare::shootArray( HepRandomEngin 49 void RandChiSquare::shootArray( HepRandomEngine* anEngine,
53 const int size, do 50 const int size, double* vect,
54 double a ) 51 double a )
55 { 52 {
56 for( double* v = vect; v != vect+size; ++v ) 53 for( double* v = vect; v != vect+size; ++v )
57 *v = shoot(anEngine,a); 54 *v = shoot(anEngine,a);
58 } 55 }
59 56
60 void RandChiSquare::fireArray( const int size, 57 void RandChiSquare::fireArray( const int size, double* vect) {
61 for( double* v = vect; v != vect+size; ++v ) 58 for( double* v = vect; v != vect+size; ++v )
62 *v = fire(defaultA); 59 *v = fire(defaultA);
63 } 60 }
64 61
65 void RandChiSquare::fireArray( const int size, 62 void RandChiSquare::fireArray( const int size, double* vect,
66 double a ) { 63 double a ) {
67 for( double* v = vect; v != vect+size; ++v ) 64 for( double* v = vect; v != vect+size; ++v )
68 *v = fire(a); 65 *v = fire(a);
69 } 66 }
70 67
71 double RandChiSquare::genChiSquare( HepRandomE 68 double RandChiSquare::genChiSquare( HepRandomEngine *anEngine,
72 double 69 double a ) {
73 /********************************************* 70 /******************************************************************
74 * 71 * *
75 * Chi Distribution - Ratio of Uniforms 72 * Chi Distribution - Ratio of Uniforms with shift *
76 * 73 * *
77 ********************************************* 74 ******************************************************************
78 * 75 * *
79 * FUNCTION : - chru samples a random number 76 * FUNCTION : - chru samples a random number from the Chi *
80 * distribution with parameter 77 * distribution with parameter a > 1. *
81 * REFERENCE : - J.F. Monahan (1987): An algo 78 * REFERENCE : - J.F. Monahan (1987): An algorithm for *
82 * generating chi random variab 79 * generating chi random variables, ACM Trans. *
83 * Math. Software 13, 168-172. 80 * Math. Software 13, 168-172. *
84 * SUBPROGRAM : - anEngine ... pointer to a ( 81 * SUBPROGRAM : - anEngine ... pointer to a (0,1)-Uniform *
85 * engine 82 * engine *
86 * 83 * *
87 * Implemented by R. Kremer, 1990 84 * Implemented by R. Kremer, 1990 *
88 ********************************************* 85 ******************************************************************/
89 86
90 static CLHEP_THREAD_LOCAL double a_in = -1.0, << 87 static double a_in = -1.0,b,vm,vp,vd;
91 double u,v,z,zz,r; 88 double u,v,z,zz,r;
92 89
93 // Check for invalid input value 90 // Check for invalid input value
94 91
95 if( a < 1 ) return (-1.0); 92 if( a < 1 ) return (-1.0);
96 93
97 if (a == 1) 94 if (a == 1)
98 { 95 {
99 for(;;) 96 for(;;)
100 { 97 {
101 u = anEngine->flat(); 98 u = anEngine->flat();
102 v = anEngine->flat() * 0.857763884960707; 99 v = anEngine->flat() * 0.857763884960707;
103 z = v / u; 100 z = v / u;
104 if (z < 0) continue; 101 if (z < 0) continue;
105 zz = z * z; 102 zz = z * z;
106 r = 2.5 - zz; 103 r = 2.5 - zz;
107 if (z < 0.0) r = r + zz * z / (3.0 * z); 104 if (z < 0.0) r = r + zz * z / (3.0 * z);
108 if (u < r * 0.3894003915) return(z*z); 105 if (u < r * 0.3894003915) return(z*z);
109 if (zz > (1.036961043 / u + 1.4)) continu 106 if (zz > (1.036961043 / u + 1.4)) continue;
110 if (2 * std::log(u) < (- zz * 0.5 )) retu 107 if (2 * std::log(u) < (- zz * 0.5 )) return(z*z);
111 } 108 }
112 } 109 }
113 else 110 else
114 { 111 {
115 if (a != a_in) 112 if (a != a_in)
116 { 113 {
117 b = std::sqrt(a - 1.0); 114 b = std::sqrt(a - 1.0);
118 vm = - 0.6065306597 * (1.0 - 0.25 / (b * 115 vm = - 0.6065306597 * (1.0 - 0.25 / (b * b + 1.0));
119 vm = (-b > vm)? -b : vm; 116 vm = (-b > vm)? -b : vm;
120 vp = 0.6065306597 * (0.7071067812 + b) / 117 vp = 0.6065306597 * (0.7071067812 + b) / (0.5 + b);
121 vd = vp - vm; 118 vd = vp - vm;
122 a_in = a; 119 a_in = a;
123 } 120 }
124 for(;;) 121 for(;;)
125 { 122 {
126 u = anEngine->flat(); 123 u = anEngine->flat();
127 v = anEngine->flat() * vd + vm; 124 v = anEngine->flat() * vd + vm;
128 z = v / u; 125 z = v / u;
129 if (z < -b) continue; 126 if (z < -b) continue;
130 zz = z * z; 127 zz = z * z;
131 r = 2.5 - zz; 128 r = 2.5 - zz;
132 if (z < 0.0) r = r + zz * z / (3.0 * (z + 129 if (z < 0.0) r = r + zz * z / (3.0 * (z + b));
133 if (u < r * 0.3894003915) return((z + b)* 130 if (u < r * 0.3894003915) return((z + b)*(z + b));
134 if (zz > (1.036961043 / u + 1.4)) continu 131 if (zz > (1.036961043 / u + 1.4)) continue;
135 if (2 * std::log(u) < (std::log(1.0 + z / 132 if (2 * std::log(u) < (std::log(1.0 + z / b) * b * b - zz * 0.5 - z * b)) return((z + b)*(z + b));
136 } 133 }
137 } 134 }
138 } 135 }
139 136
140 std::ostream & RandChiSquare::put ( std::ostre 137 std::ostream & RandChiSquare::put ( std::ostream & os ) const {
141 long pr=os.precision(20); << 138 int pr=os.precision(20);
142 std::vector<unsigned long> t(2); 139 std::vector<unsigned long> t(2);
143 os << " " << name() << "\n"; 140 os << " " << name() << "\n";
144 os << "Uvec" << "\n"; 141 os << "Uvec" << "\n";
145 t = DoubConv::dto2longs(defaultA); 142 t = DoubConv::dto2longs(defaultA);
146 os << defaultA << " " << t[0] << " " << t[1] 143 os << defaultA << " " << t[0] << " " << t[1] << "\n";
147 os.precision(pr); 144 os.precision(pr);
148 return os; 145 return os;
149 } 146 }
150 147
151 std::istream & RandChiSquare::get ( std::istre 148 std::istream & RandChiSquare::get ( std::istream & is ) {
152 std::string inName; 149 std::string inName;
153 is >> inName; 150 is >> inName;
154 if (inName != name()) { 151 if (inName != name()) {
155 is.clear(std::ios::badbit | is.rdstate()); 152 is.clear(std::ios::badbit | is.rdstate());
156 std::cerr << "Mismatch when expecting to r 153 std::cerr << "Mismatch when expecting to read state of a "
157 << name() << " distribution\n" 154 << name() << " distribution\n"
158 << "Name found was " << inName 155 << "Name found was " << inName
159 << "\nistream is left in the badbit st 156 << "\nistream is left in the badbit state\n";
160 return is; 157 return is;
161 } 158 }
162 if (possibleKeywordInput(is, "Uvec", default 159 if (possibleKeywordInput(is, "Uvec", defaultA)) {
163 std::vector<unsigned long> t(2); 160 std::vector<unsigned long> t(2);
164 is >> defaultA >> t[0] >> t[1]; defaultA = 161 is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t);
165 return is; 162 return is;
166 } 163 }
167 // is >> defaultA encompassed by possibleKey 164 // is >> defaultA encompassed by possibleKeywordInput
168 return is; 165 return is;
169 } 166 }
170 167
171 168
172 169
173 } // namespace CLHEP 170 } // namespace CLHEP
174 171
175 172