15 #include "CLHEP/Random/RandGaussQ.h"
16 #include "CLHEP/Units/PhysicalConstants.h"
25 RandGaussQ::~RandGaussQ() {
28 double RandGaussQ::operator()() {
29 return transformQuick(localEngine->flat()) * defaultStdDev + defaultMean;
32 double RandGaussQ::operator()(
double mean,
double stdDev ) {
33 return transformQuick(localEngine->flat()) * stdDev + mean;
36 void RandGaussQ::shootArray(
const int size,
double* vect,
37 double mean,
double stdDev )
39 for(
double* v = vect; v != vect + size; ++v )
40 *v =
shoot(mean,stdDev);
43 void RandGaussQ::shootArray( HepRandomEngine* anEngine,
44 const int size,
double* vect,
45 double mean,
double stdDev )
47 for(
double* v = vect; v != vect + size; ++v )
48 *v =
shoot(anEngine,mean,stdDev);
51 void RandGaussQ::fireArray(
const int size,
double* vect)
53 for(
double* v = vect; v != vect + size; ++v )
54 *v = fire( defaultMean, defaultStdDev );
57 void RandGaussQ::fireArray(
const int size,
double* vect,
58 double mean,
double stdDev )
60 for(
double* v = vect; v != vect + size; ++v )
61 *v = fire( mean, stdDev );
74 #define Table0size 250
75 #define Table1size 1000
76 #define TableSize (Table0size+Table1size)
78 #define Table0step (2.0E-6)
79 #define Table1step (5.0E-4)
81 #define Table0scale (1.0/Table1step)
83 #define Table0offset 0
84 #define Table1offset (Table0size)
89 #include "CLHEP/Random/gaussQTables.cdat"
92 double RandGaussQ::transformQuick (
double r) {
118 double y0 = gaussTables [index++];
119 double y1 = gaussTables [index];
121 return (
float) (sign * ( y1 * dx + y0 * (1.0-dx) ));
150 for (
int i = 1; i < 50; i++ ) {
151 double vn2 = 1.0/(guess*guess);
152 double s1 = -13*11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2*vn2;
153 s1 += 11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2;
154 s1 += -9*7*5*3 * vn2*vn2*vn2*vn2*vn2;
155 s1 += 7*5*3 * vn2*vn2*vn2*vn2;
156 s1 += -5*3 * vn2*vn2*vn2;
157 s1 += 3 * vn2*vn2 - vn2 + 1.0;
158 v = std::sqrt ( 2.0 * std::log ( s1 / (r*guess*std::sqrt(
CLHEP::twopi)) ) );
159 if ( std::fabs(v-guess) < eps )
break;
166 std::ostream & RandGaussQ::put ( std::ostream & os )
const {
167 int pr=os.precision(20);
168 os <<
" " <<
name() <<
"\n";
174 std::istream & RandGaussQ::get ( std::istream & is ) {
177 if (inName !=
name()) {
178 is.clear(std::ios::badbit | is.rdstate());
179 std::cerr <<
"Mismatch when expecting to read state of a "
180 <<
name() <<
" distribution\n"
181 <<
"Name found was " << inName
182 <<
"\nistream is left in the badbit state\n";
ThreeVector shoot(const G4int Ap, const G4int Af)
static int engine(pchar, pchar, double &, pchar &, const dic_type &)
static const G4double eps
static const double twopi
static const double gaussTables[2 *TableSize]
double transformSmall(double r)
G4int sign(const T t)
A simple sign function that allows us to port fortran code to c++ more easily.