CLHEP/RandomObjects/RandMultiGauss.h

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 

Class Library for High Energy Physics (version 1.8)