CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

RandStudentT.cc
Go to the documentation of this file.
1 // $Id: RandStudentT.cc,v 1.6 2010/06/16 17:24:53 garren Exp $
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()
22 #include "CLHEP/Random/defs.h"
23 #include "CLHEP/Random/RandStudentT.h"
24 #include "CLHEP/Random/DoubConv.hh"
25 
26 namespace CLHEP {
27 
28 std::string RandStudentT::name() const {return "RandStudentT";}
29 HepRandomEngine & RandStudentT::engine() {return *localEngine;}
30 
32 {
33 }
34 
36  return fire( defaultA );
37 }
38 
39 double RandStudentT::operator()( double a ) {
40  return fire( a );
41 }
42 
43 double RandStudentT::shoot( double a ) {
44 /******************************************************************
45  * *
46  * Student-t Distribution - Polar Method *
47  * *
48  ******************************************************************
49  * *
50  * The polar method of Box/Muller for generating Normal variates *
51  * is adapted to the Student-t distribution. The two generated *
52  * variates are not independent and the expected no. of uniforms *
53  * per variate is 2.5464. *
54  * *
55  ******************************************************************
56  * *
57  * FUNCTION : - tpol samples a random number from the *
58  * Student-t distribution with parameter a > 0. *
59  * REFERENCE : - R.W. Bailey (1994): Polar generation of random *
60  * variates with the t-distribution, Mathematics *
61  * of Computation 62, 779-781. *
62  * SUBPROGRAM : - ... (0,1)-Uniform generator *
63  * *
64  * *
65  * Implemented by F. Niederl, 1994 *
66  ******************************************************************/
67 
68  double u,v,w;
69 
70 // Check for valid input value
71 
72  if( a < 0.0) return (DBL_MAX);
73 
74  do
75  {
76  u = 2.0 * HepRandom::getTheEngine()->flat() - 1.0;
77  v = 2.0 * HepRandom::getTheEngine()->flat() - 1.0;
78  }
79  while ((w = u * u + v * v) > 1.0);
80 
81  return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
82 }
83 
84 void RandStudentT::shootArray( const int size, double* vect,
85  double a )
86 {
87  for( double* v = vect; v != vect + size; ++v )
88  *v = shoot(a);
89 }
90 
92  const int size, double* vect,
93  double a )
94 {
95  for( double* v = vect; v != vect + size; ++v )
96  *v = shoot(anEngine,a);
97 }
98 
99 double RandStudentT::fire( double a ) {
100 
101  double u,v,w;
102 
103  do
104  {
105  u = 2.0 * localEngine->flat() - 1.0;
106  v = 2.0 * localEngine->flat() - 1.0;
107  }
108  while ((w = u * u + v * v) > 1.0);
109 
110  return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
111 }
112 
113 void RandStudentT::fireArray( const int size, double* vect)
114 {
115  for( double* v = vect; v != vect + size; ++v )
116  *v = fire(defaultA);
117 }
118 
119 void RandStudentT::fireArray( const int size, double* vect,
120  double a )
121 {
122  for( double* v = vect; v != vect + size; ++v )
123  *v = fire(a);
124 }
125 
126 double RandStudentT::shoot( HepRandomEngine *anEngine, double a ) {
127 
128  double u,v,w;
129 
130  do
131  {
132  u = 2.0 * anEngine->flat() - 1.0;
133  v = 2.0 * anEngine->flat() - 1.0;
134  }
135  while ((w = u * u + v * v) > 1.0);
136 
137  return(u * std::sqrt( a * ( std::exp(- 2.0 / a * std::log(w)) - 1.0) / w));
138 }
139 
140 std::ostream & RandStudentT::put ( std::ostream & os ) const {
141  int 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] << "\n";
147  os.precision(pr);
148  return os;
149 #ifdef REMOVED
150  int pr=os.precision(20);
151  os << " " << name() << "\n";
152  os << defaultA << "\n";
153  os.precision(pr);
154  return os;
155 #endif
156 }
157 
158 std::istream & RandStudentT::get ( std::istream & is ) {
159  std::string inName;
160  is >> inName;
161  if (inName != name()) {
162  is.clear(std::ios::badbit | is.rdstate());
163  std::cerr << "Mismatch when expecting to read state of a "
164  << name() << " distribution\n"
165  << "Name found was " << inName
166  << "\nistream is left in the badbit state\n";
167  return is;
168  }
169  if (possibleKeywordInput(is, "Uvec", defaultA)) {
170  std::vector<unsigned long> t(2);
171  is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t);
172  return is;
173  }
174  // is >> defaultA encompassed by possibleKeywordInput
175  return is;
176 }
177 
178 } // namespace CLHEP
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
virtual double flat()=0
static HepRandomEngine * getTheEngine()
Definition: Random.cc:166
static void shootArray(const int size, double *vect, double a=1.0)
Definition: RandStudentT.cc:84
virtual ~RandStudentT()
Definition: RandStudentT.cc:31
std::string name() const
Definition: RandStudentT.cc:28
void fireArray(const int size, double *vect)
static double shoot()
HepRandomEngine & engine()
Definition: RandStudentT.cc:29
std::istream & get(std::istream &is)
std::ostream & put(std::ostream &os) const
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
@ a