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
CLHEP::RandMultiGauss::operator()
HepVector operator()()
Definition: RandMultiGauss.cc:288
CLHEP::RandMultiGauss::RandMultiGauss
RandMultiGauss(HepRandomEngine &anEngine, const HepVector &mu, const HepSymMatrix &S)
Definition: RandMultiGauss.cc:58
mu
the goal is to keep the overall false rejection probability down at the to level For each validated we discuss which of course is by necessity relative timing We take the time for a single random via one of the fastest good and at any rate the ratios will vary by around depending on the processor and memory configuration used A timing for a distribution of units would mean no time used beyond the uniform random Summary Distribution Validated Validation Rejected Past N RandGauss RandGaussT RandGaussQ bins stepwise bins RandPoisson RandPoissonT mu< 100 N=50, 000, 000 ------- mu > RandPoissonQ mu< 100 N=50, 000, 000 -------(same as RandPoissonT) mu > RandGauss shoot() etc 2.5 units Validation tests applied and a very accurate series is substituted for the table method shoot() etc 1.7 units Validation tests applied and below about **a series approximation is but we have applied this test with N up to and never get anywhere near the rejection point Analytical considerations indicate that the validation test would not reject until O(10 **24) samples were inspected. ---------------------------------------------------------- 2. RandGeneral Method since we wish to have good resolution power even if just one of the mu values is it would be unwise to dial this down any further Validation for several values o mu)
Definition: validation.doc:263
RandomVector.h
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
CLHEP
Definition: ClhepVersion.h:13
CLHEP::RandMultiGauss::~RandMultiGauss
virtual ~RandMultiGauss()
Definition: RandMultiGauss.cc:124
CLHEP::RandMultiGauss::fire
HepVector fire()
Definition: RandMultiGauss.cc:210
CLHEP::RandMultiGauss::fireArray
void fireArray(const int size, HepVector *array)
Definition: RandMultiGauss.cc:245
defs.h