| Geant4 Cross Reference |
>> 1 // $Id:$
1 // -*- C++ -*- 2 // -*- C++ -*-
2 // 3 //
3 // ------------------------------------------- 4 // -----------------------------------------------------------------------
4 // HEP Random 5 // HEP Random
5 // --- RandPoissonQ -- 6 // --- RandPoissonQ ---
6 // class header file 7 // class header file
7 // ------------------------------------------- 8 // -----------------------------------------------------------------------
8 9
9 // Class defining RandPoissonQ, which is deriv 10 // Class defining RandPoissonQ, which is derived from RandPoison.
10 // The user interface is identical; but RandGa 11 // The user interface is identical; but RandGaussQ is much faster in all cases
11 // and a bit less accurate when mu > 100. 12 // and a bit less accurate when mu > 100.
12 13
13 // =========================================== 14 // =======================================================================
14 // M. Fischler - Created: 4th Feb 2000 15 // M. Fischler - Created: 4th Feb 2000
15 // M Fischler - put and get to/from strea 16 // M Fischler - put and get to/from streams 12/10/04
16 // 17 //
17 // =========================================== 18 // =======================================================================
18 19
19 #ifndef RandPoissonQ_h 20 #ifndef RandPoissonQ_h
20 #define RandPoissonQ_h 1 21 #define RandPoissonQ_h 1
21 22
22 #include "CLHEP/Random/Random.h" 23 #include "CLHEP/Random/Random.h"
23 #include "CLHEP/Random/RandPoisson.h" 24 #include "CLHEP/Random/RandPoisson.h"
24 25
25 namespace CLHEP { 26 namespace CLHEP {
26 27
27 /** 28 /**
28 * @author 29 * @author
29 * @ingroup random 30 * @ingroup random
30 */ 31 */
31 class RandPoissonQ : public RandPoisson { 32 class RandPoissonQ : public RandPoisson {
32 33
33 public: 34 public:
34 35
35 inline RandPoissonQ ( HepRandomEngine& anEng 36 inline RandPoissonQ ( HepRandomEngine& anEngine, double b1=1.0 );
36 inline RandPoissonQ ( HepRandomEngine* anEng 37 inline RandPoissonQ ( HepRandomEngine* anEngine, double b1=1.0 );
37 // These constructors should be used to inst 38 // These constructors should be used to instantiate a RandPoissonQ
38 // distribution object defining a local engi 39 // distribution object defining a local engine for it.
39 // The static generator will be skipped usin 40 // The static generator will be skipped using the non-static methods
40 // defined below. 41 // defined below.
41 // If the engine is passed by pointer the co 42 // If the engine is passed by pointer the corresponding engine object
42 // will be deleted by the RandPoissonQ destr 43 // will be deleted by the RandPoissonQ destructor.
43 // If the engine is passed by reference the 44 // If the engine is passed by reference the corresponding engine object
44 // will not be deleted by the RandPoissonQ d 45 // will not be deleted by the RandPoissonQ destructor.
45 46
46 virtual ~RandPoissonQ(); 47 virtual ~RandPoissonQ();
47 // Destructor 48 // Destructor
48 49
49 // Save and restore to/from streams 50 // Save and restore to/from streams
50 51
51 std::ostream & put ( std::ostream & os ) con 52 std::ostream & put ( std::ostream & os ) const;
52 std::istream & get ( std::istream & is ); 53 std::istream & get ( std::istream & is );
53 54
54 // Methods to generate Poisson-distributed r 55 // Methods to generate Poisson-distributed random deviates.
55 56
56 // The method used for mu <= 100 is exact, 57 // The method used for mu <= 100 is exact, and 3-7 times faster than
57 // that used by RandPoisson. 58 // that used by RandPoisson.
58 // For mu > 100 then we use a corrected ve 59 // For mu > 100 then we use a corrected version of a
59 // (quick) Gaussian approximation. Naivel 60 // (quick) Gaussian approximation. Naively that would be:
60 // 61 //
61 // Poisson(mu) ~ floor( mu + .5 + Gaussian( 62 // Poisson(mu) ~ floor( mu + .5 + Gaussian(sqrt(mu)) )
62 // 63 //
63 // but actually, that would give a slightl 64 // but actually, that would give a slightly incorrect sigma and a
64 // very different skew than a true Poisson 65 // very different skew than a true Poisson. Instead we return
65 // 66 //
66 // Poisson(mu) ~ floor( a0*mu + a1*g + a2*g 67 // Poisson(mu) ~ floor( a0*mu + a1*g + a2*g*g ) )
67 // (with g a gaussian normal) 68 // (with g a gaussian normal)
68 // 69 //
69 // where a0, a1, a2 are chosen to give the 70 // where a0, a1, a2 are chosen to give the exctly correct mean, sigma,
70 // and skew for the Poisson distribution. 71 // and skew for the Poisson distribution.
71 72
72 // Static methods to shoot random values usi 73 // Static methods to shoot random values using the static generator
73 74
74 static long shoot( double mean=1.0 ); << 75 static long shoot( double m=1.0 );
75 76
76 static void shootArray ( const int size, lo << 77 static void shootArray ( const int size, long* vect, double m=1.0 );
77 78
78 // Static methods to shoot random values us 79 // Static methods to shoot random values using a given engine
79 // by-passing the static generator. 80 // by-passing the static generator.
80 81
81 static long shoot( HepRandomEngine* anEngin << 82 static long shoot( HepRandomEngine* anEngine, double m=1.0 );
82 83
83 static void shootArray ( HepRandomEngine* a 84 static void shootArray ( HepRandomEngine* anEngine,
84 const int size, lo << 85 const int size, long* vect, double m=1.0 );
85 86
86 // Methods using the localEngine to shoot r 87 // Methods using the localEngine to shoot random values, by-passing
87 // the static generator. 88 // the static generator.
88 89
89 long fire(); 90 long fire();
90 long fire( double mean ); << 91 long fire( double m );
91 92
92 void fireArray ( const int size, long* vect 93 void fireArray ( const int size, long* vect );
93 void fireArray ( const int size, long* vect, << 94 void fireArray ( const int size, long* vect, double m);
94 95
95 double operator()(); 96 double operator()();
96 double operator()( double mean ); << 97 double operator()( double m );
97 98
98 std::string name() const; 99 std::string name() const;
99 HepRandomEngine & engine(); 100 HepRandomEngine & engine();
100 101
101 static std::string distributionName() {retur 102 static std::string distributionName() {return "RandPoissonQ";}
102 // Provides the name of this distribution cl 103 // Provides the name of this distribution class
103 104
104 105
105 // static constants of possible interest to 106 // static constants of possible interest to users:
106 107
107 // RandPoisson will never return a deviate g 108 // RandPoisson will never return a deviate greater than this value:
108 static const double MAXIMUM_POISSON_DEVIATE; 109 static const double MAXIMUM_POISSON_DEVIATE; // Will be 2.0E9
109 110
110 static inline int tableBoundary(); 111 static inline int tableBoundary();
111 112
112 private: 113 private:
113 114
114 // constructor helper 115 // constructor helper
115 void setupForDefaultMu(); 116 void setupForDefaultMu();
116 117
117 // algorithm helper methods - all static sin 118 // algorithm helper methods - all static since the shoot methods mayneed them
118 static long poissonDeviateSmall ( HepRandomE 119 static long poissonDeviateSmall ( HepRandomEngine * e, double mean );
119 static long poissonDeviateQuick ( HepRandomE 120 static long poissonDeviateQuick ( HepRandomEngine * e, double mean );
120 static long poissonDeviateQuick ( HepRandomE 121 static long poissonDeviateQuick ( HepRandomEngine * e,
121 double A0, double A1, double A2, double si 122 double A0, double A1, double A2, double sig );
122 123
123 // All the engine info, and the default mean 124 // All the engine info, and the default mean, are in the
124 // RandPoisson base class. 125 // RandPoisson base class.
125 126
126 // quantities for approximate Poisson by cor 127 // quantities for approximate Poisson by corrected Gaussian
127 double a0; 128 double a0;
128 double a1; 129 double a1;
129 double a2; 130 double a2;
130 double sigma; 131 double sigma;
131 132
132 // static data - constants only, so that sav 133 // static data - constants only, so that saveEngineStatus works properly!
133 134
134 // The following MUST MATCH the correspondin 135 // The following MUST MATCH the corresponding values used (in
135 // poissonTables.cc) when poissonTables.cdat 136 // poissonTables.cc) when poissonTables.cdat was created.
136 // poissonTables.cc gets these values by inc 137 // poissonTables.cc gets these values by including this header,
137 // but we must be careful not to change thes 138 // but we must be careful not to change these values,
138 // and rebuild RandPoissonQ before re-genera 139 // and rebuild RandPoissonQ before re-generating poissonTables.cdat.
139 140
140 // (These statics are given values near the 141 // (These statics are given values near the start of the .cc file)
141 142
142 static const double FIRST_MU; // lowest mu 143 static const double FIRST_MU; // lowest mu value in table
143 static const double LAST_MU; // highest mu 144 static const double LAST_MU; // highest mu value
144 static const double S; // Spacing be 145 static const double S; // Spacing between mu values
145 static const int BELOW; // Starting p 146 static const int BELOW; // Starting point for N is at mu - BELOW
146 static const int ENTRIES; // Number of 147 static const int ENTRIES; // Number of entries in each mu row
147 148
148 }; 149 };
149 150
150 } // namespace CLHEP 151 } // namespace CLHEP
151 152
152 #include "CLHEP/Random/RandPoissonQ.icc" 153 #include "CLHEP/Random/RandPoissonQ.icc"
153 154
154 #endif 155 #endif
155 156