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  _mean0(right._mean0),
36  _mean1(right._mean1),
37  _mean2(right._mean2),
38  _sigma0(right._sigma0),
39  _sigma1(right._sigma1),
40  _sigma2(right._sigma2),
41  _corr01(right._corr01),
42  _corr02(right._corr02),
43  _corr12(right._corr12)
44 {
45 }
46 
47 
49  assert (a.dimension()==3);
50  double x = a[0];
51  double y = a[1];
52  double z = a[2];
53 
54 
55  double x0 = _mean0.getValue();
56  double y0 = _mean1.getValue();
57  double z0 = _mean2.getValue();
58 
59  double dx = x-x0;
60  double dy = y-y0;
61  double dz = z-z0;
62 
63  double sx = _sigma0.getValue();
64  double sy = _sigma1.getValue();
65  double sz = _sigma2.getValue();
66 
67 
68  double sxs = sx*sx;
69  double sys = sy*sy;
70  double szs = sz*sz;
71 
72 
73  double rho1 = _corr01.getValue();
74  double rho2 = _corr12.getValue();
75  double rho3 = _corr02.getValue();
76 
77  double dt = (1.0+rho1*rho2*rho3-rho1*rho1-rho2*rho2-rho3*rho3);
78  double tmp1 ,tmp2;
79 
80  tmp1= 1.0/((2*M_PI)*sqrt(2*M_PI)*sx*sy*sz*sqrt(dt));
81  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));
82 
83 
84  return tmp1*tmp2;
85 
86 
87 }
88 
90  return _mean0;
91 }
92 
94  return _sigma0;
95 }
96 
98  return _mean0;
99 }
100 
102  return _sigma0;
103 }
104 
106  return _mean1;
107 }
108 
110  return _sigma1;
111 }
112 
114  return _mean1;
115 }
116 
118  return _sigma1;
119 }
120 
122  return _mean2;
123 }
124 
125 
127  return _sigma2;
128 }
129 
131  return _mean2;
132 }
133 
135  return _sigma2;
136 }
137 
138 
139 
141  return _corr01;
142 }
143 
145  return _corr01;
146 }
147 
149  return _corr02;
150 }
151 
153  return _corr02;
154 }
155 
157  return _corr12;
158 }
159 
161  return _corr12;
162 }
163 
164 
166  return 3;
167 }
168 
169 double TrivariateGaussian::operator ()(double x) const
170 {
171  std::cerr
172  << "Warning. trivariate Gaussian called with scalar argument"
173  << std::endl;
174  assert(0);
175  return 0;
176 }
177 
178 } // namespace Genfun