18 #include "CLHEP/Random/RandChiSquare.h"
19 #include "CLHEP/Random/DoubConv.h"
20 #include "CLHEP/Utility/thread_local.h"
28 RandChiSquare::~RandChiSquare() {
32 return genChiSquare( anEngine, a );
36 HepRandomEngine *anEngine = HepRandom::getTheEngine();
37 return genChiSquare( anEngine, a );
40 double RandChiSquare::fire(
double a ) {
41 return genChiSquare( localEngine.get(),
a );
44 void RandChiSquare::shootArray(
const int size,
double* vect,
46 for(
double* v = vect; v != vect+size; ++v )
50 void RandChiSquare::shootArray( HepRandomEngine* anEngine,
51 const int size,
double* vect,
54 for(
double* v = vect; v != vect+size; ++v )
55 *v =
shoot(anEngine,a);
58 void RandChiSquare::fireArray(
const int size,
double* vect) {
59 for(
double* v = vect; v != vect+size; ++v )
63 void RandChiSquare::fireArray(
const int size,
double* vect,
65 for(
double* v = vect; v != vect+size; ++v )
69 double RandChiSquare::genChiSquare( HepRandomEngine *anEngine,
88 static CLHEP_THREAD_LOCAL
double a_in = -1.0,b,vm,vp,vd;
93 if( a < 1 )
return (-1.0);
100 v = anEngine->flat() * 0.857763884960707;
105 if (z < 0.0) r = r + zz * z / (3.0 *
z);
106 if (u < r * 0.3894003915)
return(z*z);
107 if (zz > (1.036961043 / u + 1.4))
continue;
108 if (2 * std::log(u) < (- zz * 0.5 ))
return(z*z);
115 b = std::sqrt(a - 1.0);
116 vm = - 0.6065306597 * (1.0 - 0.25 / (b * b + 1.0));
117 vm = (-b > vm)? -b : vm;
118 vp = 0.6065306597 * (0.7071067812 + b) / (0.5 + b);
124 u = anEngine->flat();
125 v = anEngine->flat() * vd + vm;
127 if (z < -b)
continue;
130 if (z < 0.0) r = r + zz * z / (3.0 * (z + b));
131 if (u < r * 0.3894003915)
return((z + b)*(z + b));
132 if (zz > (1.036961043 / u + 1.4))
continue;
133 if (2 * std::log(u) < (std::log(1.0 + z / b) * b * b - zz * 0.5 - z * b))
return((z + b)*(z + b));
138 std::ostream & RandChiSquare::put ( std::ostream & os )
const {
139 int pr=os.precision(20);
140 std::vector<unsigned long> t(2);
141 os <<
" " <<
name() <<
"\n";
142 os <<
"Uvec" <<
"\n";
143 t = DoubConv::dto2longs(defaultA);
144 os << defaultA <<
" " << t[0] <<
" " << t[1] <<
"\n";
149 std::istream & RandChiSquare::get ( std::istream & is ) {
152 if (inName !=
name()) {
153 is.clear(std::ios::badbit | is.rdstate());
154 std::cerr <<
"Mismatch when expecting to read state of a "
155 <<
name() <<
" distribution\n"
156 <<
"Name found was " << inName
157 <<
"\nistream is left in the badbit state\n";
160 if (possibleKeywordInput(is,
"Uvec", defaultA)) {
161 std::vector<unsigned long> t(2);
162 is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t);
ThreeVector shoot(const G4int Ap, const G4int Af)
static int engine(pchar, pchar, double &, pchar &, const dic_type &)