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

TrivariateGaussian.cc
Go to the documentation of this file.
1// -*- C++ -*-
2// $Id: TrivariateGaussian.cc,v 1.8 2010/06/16 18:22:01 garren Exp $
3// ---------------------------------------------------------------------------
4
7#include <assert.h>
8#include <cmath> // for exp()
9
10#if (defined __STRICT_ANSI__) || (defined _WIN32)
11#ifndef M_PI
12#define M_PI 3.14159265358979323846
13#endif // M_PI
14#endif // __STRICT_ANSI__
15
16namespace Genfun {
17FUNCTION_OBJECT_IMP(TrivariateGaussian)
18
20 _mean0("Mean0", 0.0,-10,10),
21 _mean1("Mean1", 0.0,-10,10),
22 _mean2("Mean2", 0.0,-10,10),
23 _sigma0("Sigma0",1.0,0, 10),
24 _sigma1("Sigma1",1.0,0, 10),
25 _sigma2("Sigma2",1.0,0, 10),
26 _corr01("Corr01", 0.0, -1.0, 1.0),
27 _corr02("Corr02", 0.0, -1.0, 1.0),
28 _corr12("Corr12", 0.0, -1.0, 1.0)
29{}
30
32}
33
35 AbsFunction(right),
36 _mean0(right._mean0),
37 _mean1(right._mean1),
38 _mean2(right._mean2),
39 _sigma0(right._sigma0),
40 _sigma1(right._sigma1),
41 _sigma2(right._sigma2),
42 _corr01(right._corr01),
43 _corr02(right._corr02),
44 _corr12(right._corr12)
45{
46}
47
48
50 assert (a.dimension()==3);
51 double x = a[0];
52 double y = a[1];
53 double z = a[2];
54
55
56 double x0 = _mean0.getValue();
57 double y0 = _mean1.getValue();
58 double z0 = _mean2.getValue();
59
60 double dx = x-x0;
61 double dy = y-y0;
62 double dz = z-z0;
63
64 double sx = _sigma0.getValue();
65 double sy = _sigma1.getValue();
66 double sz = _sigma2.getValue();
67
68
69 double sxs = sx*sx;
70 double sys = sy*sy;
71 double szs = sz*sz;
72
73
74 double rho1 = _corr01.getValue();
75 double rho2 = _corr12.getValue();
76 double rho3 = _corr02.getValue();
77
78 double dt = (1.0+rho1*rho2*rho3-rho1*rho1-rho2*rho2-rho3*rho3);
79 double tmp1 ,tmp2;
80
81 tmp1= 1.0/((2*M_PI)*sqrt(2*M_PI)*sx*sy*sz*sqrt(dt));
82 tmp2= exp(-0.5/dt*(dx*dx*(1.0-rho2*rho2)/sxs+dy*dy*(1.0-rho3*rho3)/sys+dz*dz*(1.0-rho1*rho1)/szs+2.0*dx*dy*(rho2*rho3-rho1)/sx/sy+2.0*dy*dz*(rho1*rho3-rho2)/sy/sz+2.0*dx*dz*(rho1*rho2-rho3)/sx/sz));
83
84
85 return tmp1*tmp2;
86
87
88}
89
91 return _mean0;
92}
93
95 return _sigma0;
96}
97
99 return _mean0;
100}
101
103 return _sigma0;
104}
105
107 return _mean1;
108}
109
111 return _sigma1;
112}
113
115 return _mean1;
116}
117
119 return _sigma1;
120}
121
123 return _mean2;
124}
125
126
128 return _sigma2;
129}
130
132 return _mean2;
133}
134
136 return _sigma2;
137}
138
139
140
142 return _corr01;
143}
144
146 return _corr01;
147}
148
150 return _corr02;
151}
152
154 return _corr02;
155}
156
158 return _corr12;
159}
160
162 return _corr12;
163}
164
165
167 return 3;
168}
169
171{
172 std::cerr
173 << "Warning. trivariate Gaussian called with scalar argument"
174 << std::endl;
175 assert(0);
176 return 0;
177}
178
179} // namespace Genfun
#define FUNCTION_OBJECT_IMP(classname)
virtual double getValue() const
Definition: Parameter.cc:27
virtual double operator()(double argument) const
virtual unsigned int dimensionality() const
@ a