CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandMultiGauss.h,v 1.3 2003/10/23 21:29:51 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandMultiGauss --- 00007 // class header file 00008 // ----------------------------------------------------------------------- 00009 00010 // Class defining methods for firing multivariate gaussian distributed 00011 // random values, given a vector of means and a covariance matrix 00012 // Definitions are those from 1998 Review of Particle Physics, section 28.3.3. 00013 // 00014 // This utilizes the following other comonents of CLHEP: 00015 // RandGauss from the Random package to get individual deviates 00016 // HepVector, HepSymMatrix and HepMatrix from the Matrix package 00017 // HepSymMatrix::similarity(HepMatrix) 00018 // diagonalize(HepSymMatrix *s) 00019 // The author of this distribution relies on diagonalize() being correct. 00020 // 00021 // Although original distribution classes in the Random package return a 00022 // double when fire() (or operator()) is done, RandMultiGauss returns a 00023 // HepVector of values. 00024 // 00025 // ======================================================================= 00026 // Mark Fischler - Created: 19th September 1999 00027 // ======================================================================= 00028 00029 #ifndef RandMultiGauss_h 00030 #define RandMultiGauss_h 1 00031 00032 #include "CLHEP/RandomObjects/defs.h" 00033 #include "CLHEP/Random/RandomEngine.h" 00034 #include "CLHEP/RandomObjects/RandomVector.h" 00035 #include "CLHEP/Matrix/Vector.h" 00036 #include "CLHEP/Matrix/Matrix.h" 00037 #include "CLHEP/Matrix/SymMatrix.h" 00038 00039 namespace CLHEP { 00040 00045 class RandMultiGauss : public HepRandomVector { 00046 00047 public: 00048 00049 RandMultiGauss ( HepRandomEngine& anEngine, 00050 const HepVector& mu, 00051 const HepSymMatrix& S ); 00052 // The symmetric matrix S MUST BE POSITIVE DEFINITE 00053 // and MUST MATCH THE SIZE OF MU. 00054 00055 RandMultiGauss ( HepRandomEngine* anEngine, 00056 const HepVector& mu, 00057 const HepSymMatrix& S ); 00058 // The symmetric matrix S MUST BE POSITIVE DEFINITE 00059 // and MUST MATCH THE SIZE OF MU. 00060 00061 // These constructors should be used to instantiate a RandMultiGauss 00062 // distribution object defining a local engine for it. 00063 // The static generator will be skipped using the non-static methods 00064 // defined below. 00065 // If the engine is passed by pointer the corresponding engine object 00066 // will be deleted by the RandMultiGauss destructor. 00067 // If the engine is passed by reference the corresponding engine object 00068 // will not be deleted by the RandGauss destructor. 00069 00070 // The following are provided for convenience in the case where each 00071 // random fired will have a different mu and S. They set the default mu 00072 // to the zero 2-vector, and the default S to the Unit 2x2 matrix. 00073 RandMultiGauss ( HepRandomEngine& anEngine ); 00074 RandMultiGauss ( HepRandomEngine* anEngine ); 00075 00076 virtual ~RandMultiGauss(); 00077 // Destructor 00078 00079 HepVector fire(); 00080 00081 HepVector fire( const HepVector& mu, const HepSymMatrix& S ); 00082 // The symmetric matrix S MUST BE POSITIVE DEFINITE 00083 // and MUST MATCH THE SIZE OF MU. 00084 00085 // A note on efficient usage when firing many vectors of Multivariate 00086 // Gaussians: For n > 2 the work needed to diagonalize S is significant. 00087 // So if you only want a collection of uncorrelated Gaussians, it will be 00088 // quicker to generate them one at a time. 00089 // 00090 // The class saves the diagonalizing matrix for the default S. 00091 // Thus generating vectors with that same S can be quite efficient. 00092 // If you require a small number of different S's, known in 00093 // advance, consider instantiating RandMulitGauss for each different S, 00094 // sharing the same engine. 00095 // 00096 // If you require a random using the default S for a distribution but a 00097 // different mu, it is most efficient to imply use the default fire() and 00098 // add the difference of the mu's to the returned HepVector. 00099 00100 void fireArray ( const int size, HepVector* array); 00101 void fireArray ( const int size, HepVector* array, 00102 const HepVector& mu, const HepSymMatrix& S ); 00103 00104 HepVector operator()(); 00105 HepVector operator()( const HepVector& mu, const HepSymMatrix& S ); 00106 // The symmetric matrix S MUST BE POSITIVE DEFINITE 00107 // and MUST MATCH THE SIZE OF MU. 00108 00109 private: 00110 00111 // Private copy constructor. Defining it here disallows use. 00112 RandMultiGauss(const RandMultiGauss &d); 00113 00114 HepRandomEngine* localEngine; 00115 bool deleteEngine; 00116 HepVector defaultMu; 00117 HepMatrix defaultU; 00118 HepVector defaultSigmas; // Each sigma is sqrt(D[i,i]) 00119 00120 bool set; 00121 double nextGaussian; 00122 00123 static void prepareUsigmas ( const HepSymMatrix & S, 00124 HepMatrix & U, 00125 HepVector & sigmas ); 00126 00127 static HepVector deviates ( const HepMatrix & U, 00128 const HepVector & sigmas, 00129 HepRandomEngine * engine, 00130 bool& available, 00131 double& next); 00132 // Returns vector of gaussian randoms based on sigmas, rotated by U, 00133 // with means of 0. 00134 00135 }; 00136 00137 } // namespace CLHEP 00138 00139 #ifdef ENABLE_BACKWARDS_COMPATIBILITY 00140 // backwards compatibility will be enabled ONLY in CLHEP 1.9 00141 using namespace CLHEP; 00142 #endif 00143 00144 #endif // RandMultiGauss_h