CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

RandomObjects/RandMultiGauss.h
Go to the documentation of this file.
1 // $Id: RandMultiGauss.h,v 1.3 2003/10/23 21:29:51 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- RandMultiGauss ---
7 // class header file
8 // -----------------------------------------------------------------------
9 
10 // Class defining methods for firing multivariate gaussian distributed
11 // random values, given a vector of means and a covariance matrix
12 // Definitions are those from 1998 Review of Particle Physics, section 28.3.3.
13 //
14 // This utilizes the following other comonents of CLHEP:
15 // RandGauss from the Random package to get individual deviates
16 // HepVector, HepSymMatrix and HepMatrix from the Matrix package
17 // HepSymMatrix::similarity(HepMatrix)
18 // diagonalize(HepSymMatrix *s)
19 // The author of this distribution relies on diagonalize() being correct.
20 //
21 // Although original distribution classes in the Random package return a
22 // double when fire() (or operator()) is done, RandMultiGauss returns a
23 // HepVector of values.
24 //
25 // =======================================================================
26 // Mark Fischler - Created: 19th September 1999
27 // =======================================================================
28 
29 #ifndef RandMultiGauss_h
30 #define RandMultiGauss_h 1
31 
33 #include "CLHEP/Random/RandomEngine.h"
35 #include "CLHEP/Matrix/Vector.h"
36 #include "CLHEP/Matrix/Matrix.h"
37 #include "CLHEP/Matrix/SymMatrix.h"
38 
39 namespace CLHEP {
40 
45 class RandMultiGauss : public HepRandomVector {
46 
47 public:
48 
49  RandMultiGauss ( HepRandomEngine& anEngine,
50  const HepVector& mu,
51  const HepSymMatrix& S );
52  // The symmetric matrix S MUST BE POSITIVE DEFINITE
53  // and MUST MATCH THE SIZE OF MU.
54 
55  RandMultiGauss ( HepRandomEngine* anEngine,
56  const HepVector& mu,
57  const HepSymMatrix& S );
58  // The symmetric matrix S MUST BE POSITIVE DEFINITE
59  // and MUST MATCH THE SIZE OF MU.
60 
61  // These constructors should be used to instantiate a RandMultiGauss
62  // distribution object defining a local engine for it.
63  // The static generator will be skipped using the non-static methods
64  // defined below.
65  // If the engine is passed by pointer the corresponding engine object
66  // will be deleted by the RandMultiGauss destructor.
67  // If the engine is passed by reference the corresponding engine object
68  // will not be deleted by the RandGauss destructor.
69 
70  // The following are provided for convenience in the case where each
71  // random fired will have a different mu and S. They set the default mu
72  // to the zero 2-vector, and the default S to the Unit 2x2 matrix.
73  RandMultiGauss ( HepRandomEngine& anEngine );
74  RandMultiGauss ( HepRandomEngine* anEngine );
75 
76  virtual ~RandMultiGauss();
77  // Destructor
78 
79  HepVector fire();
80 
81  HepVector fire( const HepVector& mu, const HepSymMatrix& S );
82  // The symmetric matrix S MUST BE POSITIVE DEFINITE
83  // and MUST MATCH THE SIZE OF MU.
84 
85  // A note on efficient usage when firing many vectors of Multivariate
86  // Gaussians: For n > 2 the work needed to diagonalize S is significant.
87  // So if you only want a collection of uncorrelated Gaussians, it will be
88  // quicker to generate them one at a time.
89  //
90  // The class saves the diagonalizing matrix for the default S.
91  // Thus generating vectors with that same S can be quite efficient.
92  // If you require a small number of different S's, known in
93  // advance, consider instantiating RandMulitGauss for each different S,
94  // sharing the same engine.
95  //
96  // If you require a random using the default S for a distribution but a
97  // different mu, it is most efficient to imply use the default fire() and
98  // add the difference of the mu's to the returned HepVector.
99 
100  void fireArray ( const int size, HepVector* array);
101  void fireArray ( const int size, HepVector* array,
102  const HepVector& mu, const HepSymMatrix& S );
103 
104  HepVector operator()();
105  HepVector operator()( const HepVector& mu, const HepSymMatrix& S );
106  // The symmetric matrix S MUST BE POSITIVE DEFINITE
107  // and MUST MATCH THE SIZE OF MU.
108 
109 private:
110 
111  // Private copy constructor. Defining it here disallows use.
112  RandMultiGauss(const RandMultiGauss &d);
113 
114  HepRandomEngine* localEngine;
115  bool deleteEngine;
116  HepVector defaultMu;
117  HepMatrix defaultU;
118  HepVector defaultSigmas; // Each sigma is sqrt(D[i,i])
119 
120  bool set;
121  double nextGaussian;
122 
123  static void prepareUsigmas ( const HepSymMatrix & S,
124  HepMatrix & U,
125  HepVector & sigmas );
126 
127  static HepVector deviates ( const HepMatrix & U,
128  const HepVector & sigmas,
129  HepRandomEngine * engine,
130  bool& available,
131  double& next);
132  // Returns vector of gaussian randoms based on sigmas, rotated by U,
133  // with means of 0.
134 
135 };
136 
137 } // namespace CLHEP
138 
139 #ifdef ENABLE_BACKWARDS_COMPATIBILITY
140 // backwards compatibility will be enabled ONLY in CLHEP 1.9
141 using namespace CLHEP;
142 #endif
143 
144 #endif // RandMultiGauss_h
RandMultiGauss(HepRandomEngine &anEngine, const HepVector &mu, const HepSymMatrix &S)
void fireArray(const int size, HepVector *array)