00001 // $Id: RandMultiGauss.h,v 1.5 2002/04/12 15:02:45 evc 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/Random/RandomEngine.h" 00033 #include "CLHEP/RandomObjects/RandomVector.h" 00034 #include "CLHEP/Matrix/Vector.h" 00035 #include "CLHEP/Matrix/Matrix.h" 00036 #include "CLHEP/Matrix/SymMatrix.h" 00037 00042 class RandMultiGauss : public HepRandomVector { 00043 00044 public: 00045 00046 RandMultiGauss ( HepRandomEngine& anEngine, 00047 const HepVector& mu, 00048 const HepSymMatrix& S ); 00049 // The symmetric matrix S MUST BE POSITIVE DEFINITE 00050 // and MUST MATCH THE SIZE OF MU. 00051 00052 RandMultiGauss ( HepRandomEngine* anEngine, 00053 const HepVector& mu, 00054 const HepSymMatrix& S ); 00055 // The symmetric matrix S MUST BE POSITIVE DEFINITE 00056 // and MUST MATCH THE SIZE OF MU. 00057 00058 // These constructors should be used to instantiate a RandMultiGauss 00059 // distribution object defining a local engine for it. 00060 // The static generator will be skipped using the non-static methods 00061 // defined below. 00062 // If the engine is passed by pointer the corresponding engine object 00063 // will be deleted by the RandMultiGauss destructor. 00064 // If the engine is passed by reference the corresponding engine object 00065 // will not be deleted by the RandGauss destructor. 00066 00067 // The following are provided for convenience in the case where each 00068 // random fired will have a different mu and S. They set the default mu 00069 // to the zero 2-vector, and the default S to the Unit 2x2 matrix. 00070 RandMultiGauss ( HepRandomEngine& anEngine ); 00071 RandMultiGauss ( HepRandomEngine* anEngine ); 00072 00073 virtual ~RandMultiGauss(); 00074 // Destructor 00075 00076 HepVector fire(); 00077 00078 HepVector fire( const HepVector& mu, const HepSymMatrix& S ); 00079 // The symmetric matrix S MUST BE POSITIVE DEFINITE 00080 // and MUST MATCH THE SIZE OF MU. 00081 00082 // A note on efficient usage when firing many vectors of Multivariate 00083 // Gaussians: For n > 2 the work needed to diagonalize S is significant. 00084 // So if you only want a collection of uncorrelated Gaussians, it will be 00085 // quicker to generate them one at a time. 00086 // 00087 // The class saves the diagonalizing matrix for the default S. 00088 // Thus generating vectors with that same S can be quite efficient. 00089 // If you require a small number of different S's, known in 00090 // advance, consider instantiating RandMulitGauss for each different S, 00091 // sharing the same engine. 00092 // 00093 // If you require a random using the default S for a distribution but a 00094 // different mu, it is most efficient to imply use the default fire() and 00095 // add the difference of the mu's to the returned HepVector. 00096 00097 void fireArray ( const int size, HepVector* array); 00098 void fireArray ( const int size, HepVector* array, 00099 const HepVector& mu, const HepSymMatrix& S ); 00100 00101 HepVector operator()(); 00102 HepVector operator()( const HepVector& mu, const HepSymMatrix& S ); 00103 // The symmetric matrix S MUST BE POSITIVE DEFINITE 00104 // and MUST MATCH THE SIZE OF MU. 00105 00106 private: 00107 00108 // Private copy constructor. Defining it here disallows use. 00109 RandMultiGauss(const RandMultiGauss &d); 00110 00111 HepRandomEngine* localEngine; 00112 bool deleteEngine; 00113 HepVector defaultMu; 00114 HepMatrix defaultU; 00115 HepVector defaultSigmas; // Each sigma is sqrt(D[i,i]) 00116 00117 bool set; 00118 double nextGaussian; 00119 00120 static void prepareUsigmas ( const HepSymMatrix & S, 00121 HepMatrix & U, 00122 HepVector & sigmas ); 00123 00124 static HepVector deviates ( const HepMatrix & U, 00125 const HepVector & sigmas, 00126 HepRandomEngine * engine, 00127 bool& available, 00128 double& next); 00129 // Returns vector of gaussian randoms based on sigmas, rotated by U, 00130 // with means of 0. 00131 00132 }; 00133 00134 #endif // RandMultiGauss_h