Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
RandStudentT.cc
Go to the documentation of this file.
1 // $Id:$
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- RandStudentT ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 
10 // =======================================================================
11 // John Marraffino - Created: 12th May 1998
12 // G.Cosmo - Fixed minor bug on inline definition for shoot()
13 // methods : 20th Aug 1998
14 // M Fischler - put and get to/from streams 12/13/04
15 // M Fischler - put/get to/from streams uses pairs of ulongs when
16 // + storing doubles avoid problems with precision
17 // 4/14/05
18 // =======================================================================
19 
20 #include <float.h>
21 #include <cmath> // for std::log() std::exp()
23 #include "CLHEP/Random/DoubConv.h"
24 
25 namespace CLHEP {
26 
27 std::string RandStudentT::name() const {return "RandStudentT";}
28 HepRandomEngine & RandStudentT::engine() {return *localEngine;}
29 
31 {
32 }
33 
35  return fire( defaultA );
36 }
37 
38 double RandStudentT::operator()( double a ) {
39  return fire( a );
40 }
41 
42 double RandStudentT::shoot( double a ) {
43 /******************************************************************
44  * *
45  * Student-t Distribution - Polar Method *
46  * *
47  ******************************************************************
48  * *
49  * The polar method of Box/Muller for generating Normal variates *
50  * is adapted to the Student-t distribution. The two generated *
51  * variates are not independent and the expected no. of uniforms *
52  * per variate is 2.5464. *
53  * *
54  ******************************************************************
55  * *
56  * FUNCTION : - tpol samples a random number from the *
57  * Student-t distribution with parameter a > 0. *
58  * REFERENCE : - R.W. Bailey (1994): Polar generation of random *
59  * variates with the t-distribution, Mathematics *
60  * of Computation 62, 779-781. *
61  * SUBPROGRAM : - ... (0,1)-Uniform generator *
62  * *
63  * *
64  * Implemented by F. Niederl, 1994 *
65  ******************************************************************/
66 
67  double u,v,w;
68 
69 // Check for valid input value
70 
71  if( a < 0.0) return (DBL_MAX);
72 
73  do
74  {
75  u = 2.0 * HepRandom::getTheEngine()->flat() - 1.0;
76  v = 2.0 * HepRandom::getTheEngine()->flat() - 1.0;
77  }
78  while ((w = u * u + v * v) > 1.0);
79 
80  return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
81 }
82 
83 void RandStudentT::shootArray( const int size, double* vect,
84  double a )
85 {
86  for( double* v = vect; v != vect + size; ++v )
87  *v = shoot(a);
88 }
89 
91  const int size, double* vect,
92  double a )
93 {
94  for( double* v = vect; v != vect + size; ++v )
95  *v = shoot(anEngine,a);
96 }
97 
98 double RandStudentT::fire( double a ) {
99 
100  double u,v,w;
101 
102  do
103  {
104  u = 2.0 * localEngine->flat() - 1.0;
105  v = 2.0 * localEngine->flat() - 1.0;
106  }
107  while ((w = u * u + v * v) > 1.0);
108 
109  return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
110 }
111 
112 void RandStudentT::fireArray( const int size, double* vect)
113 {
114  for( double* v = vect; v != vect + size; ++v )
115  *v = fire(defaultA);
116 }
117 
118 void RandStudentT::fireArray( const int size, double* vect,
119  double a )
120 {
121  for( double* v = vect; v != vect + size; ++v )
122  *v = fire(a);
123 }
124 
125 double RandStudentT::shoot( HepRandomEngine *anEngine, double a ) {
126 
127  double u,v,w;
128 
129  do
130  {
131  u = 2.0 * anEngine->flat() - 1.0;
132  v = 2.0 * anEngine->flat() - 1.0;
133  }
134  while ((w = u * u + v * v) > 1.0);
135 
136  return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
137 }
138 
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";
146  os.precision(pr);
147  return os;
148 }
149 
150 std::istream & RandStudentT::get ( std::istream & is ) {
151  std::string inName;
152  is >> inName;
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";
159  return is;
160  }
161  if (possibleKeywordInput(is, "Uvec", defaultA)) {
162  std::vector<unsigned long> t(2);
163  is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t);
164  return is;
165  }
166  // is >> defaultA encompassed by possibleKeywordInput
167  return is;
168 }
169 
170 } // namespace CLHEP
void fireArray(const int size, double *vect)
std::istream & get(std::istream &is)
std::string name() const
Definition: RandStudentT.cc:27
virtual double flat()=0
static HepRandomEngine * getTheEngine()
Definition: Random.cc:265
virtual ~RandStudentT()
Definition: RandStudentT.cc:30
static double shoot()
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
Definition: RandomEngine.h:167
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:98
static void shootArray(const int size, double *vect, double a=1.0)
Definition: RandStudentT.cc:83
std::ostream & put(std::ostream &os) const
HepRandomEngine & engine()
Definition: RandStudentT.cc:28
#define DBL_MAX
Definition: templates.hh:83
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:114