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