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 
16 namespace Genfun {
17 FUNCTION_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
virtual unsigned int dimensionality() const
virtual double operator()(double argument) const
unsigned int dimension() const
#define FUNCTION_OBJECT_IMP(classname)
virtual double getValue() const
Definition: Parameter.cc:27