00001 // $Id: RandPoissonQ.h,v 1.5 2002/04/12 15:02:45 evc Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandPoissonQ --- 00007 // class header file 00008 // ----------------------------------------------------------------------- 00009 00010 // Class defining RandPoissonQ, which is derived from RandPoison. 00011 // The user interface is identical; but RandGaussQ is much faster in all cases 00012 // and a bit less accurate when mu > 100. 00013 00014 // ======================================================================= 00015 // M. Fischler - Created: 4th Feb 2000 00016 // 00017 // ======================================================================= 00018 00019 #ifndef RandPoissonQ_h 00020 #define RandPoissonQ_h 1 00021 00022 #include "CLHEP/Random/Random.h" 00023 #include "CLHEP/Random/RandPoisson.h" 00024 00029 class RandPoissonQ : public RandPoisson { 00030 00031 public: 00032 00033 inline RandPoissonQ ( HepRandomEngine& anEngine, double m=1.0 ); 00034 inline RandPoissonQ ( HepRandomEngine* anEngine, double m=1.0 ); 00035 // These constructors should be used to instantiate a RandPoissonQ 00036 // distribution object defining a local engine for it. 00037 // The static generator will be skipped using the non-static methods 00038 // defined below. 00039 // If the engine is passed by pointer the corresponding engine object 00040 // will be deleted by the RandPoissonQ destructor. 00041 // If the engine is passed by reference the corresponding engine object 00042 // will not be deleted by the RandPoissonQ destructor. 00043 00044 virtual ~RandPoissonQ(); 00045 // Destructor 00046 00047 // Methods to generate Poisson-distributed random deviates. 00048 00049 // The method used for mu <= 100 is exact, and 3-7 times faster than 00050 // that used by RandPoisson. 00051 // For mu > 100 then we use a corrected version of a 00052 // (quick) Gaussian approximation. Naively that would be: 00053 // 00054 // Poisson(mu) ~ floor( mu + .5 + Gaussian(sqrt(mu)) ) 00055 // 00056 // but actually, that would give a slightly incorrect sigma and a 00057 // very different skew than a true Poisson. Instead we return 00058 // 00059 // Poisson(mu) ~ floor( a0*mu + a1*g + a2*g*g ) ) 00060 // (with g a gaussian normal) 00061 // 00062 // where a0, a1, a2 are chosen to give the exctly correct mean, sigma, 00063 // and skew for the Poisson distribution. 00064 00065 // Static methods to shoot random values using the static generator 00066 00067 static long shoot( double m=1.0 ); 00068 00069 static void shootArray ( const int size, long* vect, double m=1.0 ); 00070 00071 // Static methods to shoot random values using a given engine 00072 // by-passing the static generator. 00073 00074 static long shoot( HepRandomEngine* anEngine, double m=1.0 ); 00075 00076 static void shootArray ( HepRandomEngine* anEngine, 00077 const int size, long* vect, double m=1.0 ); 00078 00079 // Methods using the localEngine to shoot random values, by-passing 00080 // the static generator. 00081 00082 long fire(); 00083 long fire( double m ); 00084 00085 void fireArray ( const int size, long* vect ); 00086 void fireArray ( const int size, long* vect, double m); 00087 00088 double operator()(); 00089 double operator()( double m ); 00090 00091 // static constants of possible interest to users: 00092 00093 // RandPoisson will never return a deviate greater than this value: 00094 static const double MAXIMUM_POISSON_DEVIATE; // Will be 2.0E9 00095 00096 static inline int tableBoundary(); 00097 00098 protected: 00099 00100 private: 00101 00102 // Private copy constructor. Defining it here disallows use. 00103 RandPoissonQ(const RandPoissonQ& d); 00104 00105 // constructor helper 00106 void setupForDefaultMu(); 00107 00108 // algorithm helper methods - all static since the shoot methods mayneed them 00109 static long poissonDeviateSmall ( HepRandomEngine * e, double mean ); 00110 static long poissonDeviateQuick ( HepRandomEngine * e, double mean ); 00111 static long poissonDeviateQuick ( HepRandomEngine * e, 00112 double A0, double A1, double A2, double sig ); 00113 00114 // All the engine info, and the default mean, are in the 00115 // RandPoisson base class. 00116 00117 // quantities for approximate Poisson by corrected Gaussian 00118 double a0; 00119 double a1; 00120 double a2; 00121 double sigma; 00122 00123 // static data - constants only, so that saveEngineStatus works properly! 00124 00125 // The following MUST MATCH the corresponding values used (in 00126 // poissonTables.cc) when poissonTables.cdat was created. 00127 // poissonTables.cc gets these values by including this header, 00128 // but we must be careful not to change these values, 00129 // and rebuild RandPoissonQ before re-generating poissonTables.cdat. 00130 00131 // (These statics are given values near the start of the .cc file) 00132 00133 static const double FIRST_MU; // lowest mu value in table 00134 static const double LAST_MU; // highest mu value 00135 static const double S; // Spacing between mu values 00136 static const int BELOW; // Starting point for N is at mu - BELOW 00137 static const int ENTRIES; // Number of entries in each mu row 00138 00139 }; 00140 00141 #include "CLHEP/Random/RandPoissonQ.icc" 00142 00143 #endif