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