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

BivariateGaussian.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id: BivariateGaussian.cc,v 1.8 2010/06/16 18:22:01 garren Exp $
5 #include <assert.h>
6 #include <cmath> // for exp()
7 
8 #if (defined __STRICT_ANSI__) || (defined _WIN32)
9 #ifndef M_PI
10 #define M_PI 3.14159265358979323846
11 #endif // M_PI
12 #endif // __STRICT_ANSI__
13 
14 namespace Genfun {
15 FUNCTION_OBJECT_IMP(BivariateGaussian)
16 
18  _mean0("Mean0", 0.0,-10,10),
19  _mean1("Mean1", 0.0,-10,10),
20  _sigma0("Sigma0",1.0,0, 10),
21  _sigma1("Sigma1",1.0,0, 10),
22  _corr01("Corr01", 0.0, -1.0, 1.0)
23 {}
24 
26 }
27 
29  _mean0(right._mean0),
30  _mean1(right._mean1),
31  _sigma0(right._sigma0),
32  _sigma1(right._sigma1),
33  _corr01(right._corr01)
34 {
35 }
36 
37 double BivariateGaussian::operator() (const Argument & a) const {
38  assert (a.dimension()==2);
39  double x = a[0];
40  double y = a[1];
41 
42  double x0 = _mean0.getValue();
43  double y0 = _mean1.getValue();
44  double dx = x-x0;
45  double dy = y-y0;
46 
47  double sx = _sigma0.getValue();
48  double sy = _sigma1.getValue();
49 
50  double sxs = sx*sx;
51  double sys = sy*sy;
52  double rho = _corr01.getValue();
53  double dt = (1.0+rho)*(1.0-rho);
54 
55  return (1.0/(2*M_PI*sx*sy*sqrt(dt))) *
56  exp(-1.0/(2.0*dt)*(dx*dx/sxs+dy*dy/sys-2.0*rho*dx*dy/sx/sy));
57 }
58 
60  return _mean0;
61 }
62 
64  return _sigma0;
65 }
66 
68  return _mean0;
69 }
70 
72  return _sigma0;
73 }
74 
76  return _mean1;
77 }
78 
80  return _sigma1;
81 }
82 
84  return _mean1;
85 }
86 
88  return _sigma1;
89 }
90 
91 
92 
94  return _corr01;
95 }
96 
98  return _corr01;
99 }
100 
101 
102 unsigned int BivariateGaussian::dimensionality() const {
103  return 2;
104 }
105 
106 double BivariateGaussian::operator ()(double x) const
107 {
108  std::cerr
109  << "Warning. bivariate Gaussian called with scalar argument"
110  << std::endl;
111  assert(0);
112  return 0;
113 }
114 
115 } // namespace Genfun