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

RandChiSquare.cc
Go to the documentation of this file.
1 // $Id: RandChiSquare.cc,v 1.6 2010/06/16 17:24:53 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- RandChiSquare ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 
10 // =======================================================================
11 // John Marraffino - Created: 12th May 1998
12 // M Fischler - put and get to/from streams 12/10/04
13 // M Fischler - put/get to/from streams uses pairs of ulongs when
14 // + storing doubles avoid problems with precision
15 // 4/14/05
16 // =======================================================================
17 
18 #include "CLHEP/Random/defs.h"
19 #include "CLHEP/Random/RandChiSquare.h"
20 #include "CLHEP/Random/DoubConv.hh"
21 #include <cmath> // for std::log()
22 
23 namespace CLHEP {
24 
25 std::string RandChiSquare::name() const {return "RandChiSquare";}
26 HepRandomEngine & RandChiSquare::engine() {return *localEngine;}
27 
29 }
30 
31 double RandChiSquare::shoot( HepRandomEngine *anEngine, double a ) {
32  return genChiSquare( anEngine, a );
33 }
34 
35 double RandChiSquare::shoot( double a ) {
37  return genChiSquare( anEngine, a );
38 }
39 
40 double RandChiSquare::fire( double a ) {
41  return genChiSquare( localEngine.get(), a );
42 }
43 
44 void RandChiSquare::shootArray( const int size, double* vect,
45  double a ) {
46  for( double* v = vect; v != vect+size; ++v )
47  *v = shoot(a);
48 }
49 
51  const int size, double* vect,
52  double a )
53 {
54  for( double* v = vect; v != vect+size; ++v )
55  *v = shoot(anEngine,a);
56 }
57 
58 void RandChiSquare::fireArray( const int size, double* vect) {
59  for( double* v = vect; v != vect+size; ++v )
60  *v = fire(defaultA);
61 }
62 
63 void RandChiSquare::fireArray( const int size, double* vect,
64  double a ) {
65  for( double* v = vect; v != vect+size; ++v )
66  *v = fire(a);
67 }
68 
69 double RandChiSquare::genChiSquare( HepRandomEngine *anEngine,
70  double a ) {
71 /******************************************************************
72  * *
73  * Chi Distribution - Ratio of Uniforms with shift *
74  * *
75  ******************************************************************
76  * *
77  * FUNCTION : - chru samples a random number from the Chi *
78  * distribution with parameter a > 1. *
79  * REFERENCE : - J.F. Monahan (1987): An algorithm for *
80  * generating chi random variables, ACM Trans. *
81  * Math. Software 13, 168-172. *
82  * SUBPROGRAM : - anEngine ... pointer to a (0,1)-Uniform *
83  * engine *
84  * *
85  * Implemented by R. Kremer, 1990 *
86  ******************************************************************/
87 
88  static double a_in = -1.0,b,vm,vp,vd;
89  double u,v,z,zz,r;
90 
91 // Check for invalid input value
92 
93  if( a < 1 ) return (-1.0);
94 
95  if (a == 1)
96  {
97  for(;;)
98  {
99  u = anEngine->flat();
100  v = anEngine->flat() * 0.857763884960707;
101  z = v / u;
102  if (z < 0) continue;
103  zz = z * z;
104  r = 2.5 - zz;
105  if (z < 0.0) r = r + zz * z / (3.0 * z);
106  if (u < r * 0.3894003915) return(z*z);
107  if (zz > (1.036961043 / u + 1.4)) continue;
108  if (2 * std::log(u) < (- zz * 0.5 )) return(z*z);
109  }
110  }
111  else
112  {
113  if (a != a_in)
114  {
115  b = std::sqrt(a - 1.0);
116  vm = - 0.6065306597 * (1.0 - 0.25 / (b * b + 1.0));
117  vm = (-b > vm)? -b : vm;
118  vp = 0.6065306597 * (0.7071067812 + b) / (0.5 + b);
119  vd = vp - vm;
120  a_in = a;
121  }
122  for(;;)
123  {
124  u = anEngine->flat();
125  v = anEngine->flat() * vd + vm;
126  z = v / u;
127  if (z < -b) continue;
128  zz = z * z;
129  r = 2.5 - zz;
130  if (z < 0.0) r = r + zz * z / (3.0 * (z + b));
131  if (u < r * 0.3894003915) return((z + b)*(z + b));
132  if (zz > (1.036961043 / u + 1.4)) continue;
133  if (2 * std::log(u) < (std::log(1.0 + z / b) * b * b - zz * 0.5 - z * b)) return((z + b)*(z + b));
134  }
135  }
136 }
137 
138 std::ostream & RandChiSquare::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";
143  t = DoubConv::dto2longs(defaultA);
144  os << defaultA << " " << t[0] << " " << t[1] << "\n";
145  os.precision(pr);
146  return os;
147 #ifdef REMOVED
148  int pr=os.precision(20);
149  os << " " << name() << "\n";
150  os << defaultA << "\n";
151  os.precision(pr);
152  return os;
153 #endif
154 }
155 
156 std::istream & RandChiSquare::get ( std::istream & is ) {
157  std::string inName;
158  is >> inName;
159  if (inName != name()) {
160  is.clear(std::ios::badbit | is.rdstate());
161  std::cerr << "Mismatch when expecting to read state of a "
162  << name() << " distribution\n"
163  << "Name found was " << inName
164  << "\nistream is left in the badbit state\n";
165  return is;
166  }
167  if (possibleKeywordInput(is, "Uvec", defaultA)) {
168  std::vector<unsigned long> t(2);
169  is >> defaultA >> t[0] >> t[1]; defaultA = DoubConv::longs2double(t);
170  return is;
171  }
172  // is >> defaultA encompassed by possibleKeywordInput
173  return is;
174 }
175 
176 
177 
178 } // namespace CLHEP
179 
CLHEP::RandChiSquare::fire
double fire()
CLHEP::RandChiSquare::name
std::string name() const
Definition: RandChiSquare.cc:25
a
@ a
Definition: testCategories.cc:125
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
CLHEP::RandChiSquare::get
std::istream & get(std::istream &is)
Definition: RandChiSquare.cc:156
b
@ b
Definition: testCategories.cc:125
is
HepRotation and so forth isNear() norm2() rectify() static Rotation row1 row4(To avoid bloat in the code pulled in for programs which don 't use all these features, we split the implementation .cc files. Only isNear() goes into the original Rotation.cc) --------------------------------------- HepAxisAngle and HepEulerAngles classes --------------------------------------- These classes are very useful and simple structures for holding the result of a nice intuituve decomposition of a rotation there is no longer much content in the distinct ZOOM PhysicsVectors library The only content left in the library is the object files representing the various Exception objects When we build the CLHEP classes for the ZOOM we will set up so as to use ZOOM SpaceVector is(but we can disable namespace usage and most of our users do so at this point). What I do is leave Hep3Vector in the global namespace
CLHEP::DoubConv::longs2double
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
CLHEP::RandChiSquare::shootArray
static void shootArray(const int size, double *vect, double a=1.0)
Definition: RandChiSquare.cc:44
CLHEP::RandChiSquare::~RandChiSquare
virtual ~RandChiSquare()
Definition: RandChiSquare.cc:28
CLHEP::RandChiSquare::engine
HepRandomEngine & engine()
Definition: RandChiSquare.cc:26
CLHEP
Definition: ClhepVersion.h:13
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
CLHEP::RandChiSquare::put
std::ostream & put(std::ostream &os) const
Definition: RandChiSquare.cc:138
CLHEP::RandChiSquare::shoot
static double shoot()
CLHEP::possibleKeywordInput
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
Definition: Matrix/CLHEP/Random/RandomEngine.h:168
CLHEP::HepRandomEngine::flat
virtual double flat()=0
CLHEP::DoubConv::dto2longs
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
CLHEP::HepRandom::getTheEngine
static HepRandomEngine * getTheEngine()
Definition: Random.cc:166
CLHEP::RandChiSquare::fireArray
void fireArray(const int size, double *vect)
Definition: RandChiSquare.cc:58