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"
25
26namespace CLHEP {
27
28std::string RandStudentT::name() const {return "RandStudentT";}
29HepRandomEngine & RandStudentT::engine() {return *localEngine;}
30
34
36 return fire( defaultA );
37}
38
39double RandStudentT::operator()( double a ) {
40 return fire( a );
41}
42
43double 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
84void 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
99double 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
113void RandStudentT::fireArray( const int size, double* vect)
114{
115 for( double* v = vect; v != vect + size; ++v )
116 *v = fire(defaultA);
117}
118
119void 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
126double 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
140std::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
158std::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)
std::string name() const
void fireArray(const int size, double *vect)
static double shoot()
HepRandomEngine & engine()
std::istream & get(std::istream &is)
std::ostream & put(std::ostream &os) const
bool possibleKeywordInput(IS &is, const std::string &key, T &t)