ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RandChiSquare.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RandChiSquare.cc
1 // -*- C++ -*-
2 //
3 // -----------------------------------------------------------------------
4 // HEP Random
5 // --- RandChiSquare ---
6 // class implementation file
7 // -----------------------------------------------------------------------
8 
9 // =======================================================================
10 // John Marraffino - Created: 12th May 1998
11 // M Fischler - put and get to/from streams 12/10/04
12 // M Fischler - put/get to/from streams uses pairs of ulongs when
13 // + storing doubles avoid problems with precision
14 // 4/14/05
15 // =======================================================================
16 
18 #include "CLHEP/Random/DoubConv.h"
20 #include <cmath> // for std::log()
21 
22 namespace CLHEP {
23 
24 std::string RandChiSquare::name() const {return "RandChiSquare";}
26 
28 }
29 
30 double RandChiSquare::shoot( HepRandomEngine *anEngine, double a ) {
31  return genChiSquare( anEngine, a );
32 }
33 
34 double RandChiSquare::shoot( double a ) {
36  return genChiSquare( anEngine, a );
37 }
38 
39 double RandChiSquare::fire( double a ) {
40  return genChiSquare( localEngine.get(), a );
41 }
42 
43 void RandChiSquare::shootArray( const int size, double* vect,
44  double a ) {
45  for( double* v = vect; v != vect+size; ++v )
46  *v = shoot(a);
47 }
48 
50  const int size, double* vect,
51  double a )
52 {
53  for( double* v = vect; v != vect+size; ++v )
54  *v = shoot(anEngine,a);
55 }
56 
57 void RandChiSquare::fireArray( const int size, double* vect) {
58  for( double* v = vect; v != vect+size; ++v )
59  *v = fire(defaultA);
60 }
61 
62 void RandChiSquare::fireArray( const int size, double* vect,
63  double a ) {
64  for( double* v = vect; v != vect+size; ++v )
65  *v = fire(a);
66 }
67 
69  double a ) {
70 /******************************************************************
71  * *
72  * Chi Distribution - Ratio of Uniforms with shift *
73  * *
74  ******************************************************************
75  * *
76  * FUNCTION : - chru samples a random number from the Chi *
77  * distribution with parameter a > 1. *
78  * REFERENCE : - J.F. Monahan (1987): An algorithm for *
79  * generating chi random variables, ACM Trans. *
80  * Math. Software 13, 168-172. *
81  * SUBPROGRAM : - anEngine ... pointer to a (0,1)-Uniform *
82  * engine *
83  * *
84  * Implemented by R. Kremer, 1990 *
85  ******************************************************************/
86 
87  static CLHEP_THREAD_LOCAL double a_in = -1.0,b,vm,vp,vd;
88  double u,v,z,zz,r;
89 
90 // Check for invalid input value
91 
92  if( a < 1 ) return (-1.0);
93 
94  if (a == 1)
95  {
96  for(;;)
97  {
98  u = anEngine->flat();
99  v = anEngine->flat() * 0.857763884960707;
100  z = v / u;
101  if (z < 0) continue;
102  zz = z * z;
103  r = 2.5 - zz;
104  if (z < 0.0) r = r + zz * z / (3.0 * z);
105  if (u < r * 0.3894003915) return(z*z);
106  if (zz > (1.036961043 / u + 1.4)) continue;
107  if (2 * std::log(u) < (- zz * 0.5 )) return(z*z);
108  }
109  }
110  else
111  {
112  if (a != a_in)
113  {
114  b = std::sqrt(a - 1.0);
115  vm = - 0.6065306597 * (1.0 - 0.25 / (b * b + 1.0));
116  vm = (-b > vm)? -b : vm;
117  vp = 0.6065306597 * (0.7071067812 + b) / (0.5 + b);
118  vd = vp - vm;
119  a_in = a;
120  }
121  for(;;)
122  {
123  u = anEngine->flat();
124  v = anEngine->flat() * vd + vm;
125  z = v / u;
126  if (z < -b) continue;
127  zz = z * z;
128  r = 2.5 - zz;
129  if (z < 0.0) r = r + zz * z / (3.0 * (z + b));
130  if (u < r * 0.3894003915) return((z + b)*(z + b));
131  if (zz > (1.036961043 / u + 1.4)) continue;
132  if (2 * std::log(u) < (std::log(1.0 + z / b) * b * b - zz * 0.5 - z * b)) return((z + b)*(z + b));
133  }
134  }
135 }
136 
137 std::ostream & RandChiSquare::put ( std::ostream & os ) const {
138  int pr=os.precision(20);
139  std::vector<unsigned long> t(2);
140  os << " " << name() << "\n";
141  os << "Uvec" << "\n";
143  os << defaultA << " " << t[0] << " " << t[1] << "\n";
144  os.precision(pr);
145  return os;
146 }
147 
148 std::istream & RandChiSquare::get ( std::istream & is ) {
149  std::string inName;
150  is >> inName;
151  if (inName != name()) {
152  is.clear(std::ios::badbit | is.rdstate());
153  std::cerr << "Mismatch when expecting to read state of a "
154  << name() << " distribution\n"
155  << "Name found was " << inName
156  << "\nistream is left in the badbit state\n";
157  return is;
158  }
159  if (possibleKeywordInput(is, "Uvec", defaultA)) {
160  std::vector<unsigned long> t(2);
161  is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t);
162  return is;
163  }
164  // is >> defaultA encompassed by possibleKeywordInput
165  return is;
166 }
167 
168 
169 
170 } // namespace CLHEP
171