22 #include "CLHEP/Random/RandStudentT.h"
23 #include "CLHEP/Random/DoubConv.h"
30 RandStudentT::~RandStudentT()
34 double RandStudentT::operator()() {
35 return fire( defaultA );
38 double RandStudentT::operator()(
double a ) {
75 u = 2.0 * HepRandom::getTheEngine()->flat() - 1.0;
76 v = 2.0 * HepRandom::getTheEngine()->flat() - 1.0;
78 while ((w = u * u + v * v) > 1.0);
80 return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
83 void RandStudentT::shootArray(
const int size,
double* vect,
86 for(
double* v = vect; v != vect + size; ++v )
90 void RandStudentT::shootArray( HepRandomEngine* anEngine,
91 const int size,
double* vect,
94 for(
double* v = vect; v != vect + size; ++v )
95 *v =
shoot(anEngine,a);
98 double RandStudentT::fire(
double a ) {
104 u = 2.0 * localEngine->flat() - 1.0;
105 v = 2.0 * localEngine->flat() - 1.0;
107 while ((w = u * u + v * v) > 1.0);
109 return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
112 void RandStudentT::fireArray(
const int size,
double* vect)
114 for(
double* v = vect; v != vect + size; ++v )
118 void RandStudentT::fireArray(
const int size,
double* vect,
121 for(
double* v = vect; v != vect + size; ++v )
131 u = 2.0 * anEngine->flat() - 1.0;
132 v = 2.0 * anEngine->flat() - 1.0;
134 while ((w = u * u + v * v) > 1.0);
136 return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
139 std::ostream & RandStudentT::put ( std::ostream & os )
const {
140 int pr=os.precision(20);
141 std::vector<unsigned long> t(2);
142 os <<
" " <<
name() <<
"\n";
143 os <<
"Uvec" <<
"\n";
144 t = DoubConv::dto2longs(defaultA);
145 os << defaultA <<
" " << t[0] <<
" " << t[1] <<
"\n";
150 std::istream & RandStudentT::get ( std::istream & is ) {
153 if (inName !=
name()) {
154 is.clear(std::ios::badbit | is.rdstate());
155 std::cerr <<
"Mismatch when expecting to read state of a "
156 <<
name() <<
" distribution\n"
157 <<
"Name found was " << inName
158 <<
"\nistream is left in the badbit state\n";
161 if (possibleKeywordInput(is,
"Uvec", defaultA)) {
162 std::vector<unsigned long> t(2);
163 is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t);
ThreeVector shoot(const G4int Ap, const G4int Af)
const G4double w[NPOINTSGL]
static int engine(pchar, pchar, double &, pchar &, const dic_type &)