CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: TrivariateGaussian.cc,v 1.5.4.1.2.1 2009/11/10 20:43:12 garren Exp $ 00003 // --------------------------------------------------------------------------- 00004 00005 #include "CLHEP/GenericFunctions/defs.h" 00006 #include "CLHEP/GenericFunctions/TrivariateGaussian.hh" 00007 #include <assert.h> 00008 #include <cmath> // for exp() 00009 00010 #if (defined __STRICT_ANSI__) || (defined _WIN32) 00011 #ifndef M_PI 00012 #define M_PI 3.14159265358979323846 00013 #endif // M_PI 00014 #endif // __STRICT_ANSI__ 00015 00016 namespace Genfun { 00017 FUNCTION_OBJECT_IMP(TrivariateGaussian) 00018 00019 TrivariateGaussian::TrivariateGaussian(): 00020 _mean0("Mean0", 0.0,-10,10), 00021 _mean1("Mean1", 0.0,-10,10), 00022 _mean2("Mean2", 0.0,-10,10), 00023 _sigma0("Sigma0",1.0,0, 10), 00024 _sigma1("Sigma1",1.0,0, 10), 00025 _sigma2("Sigma2",1.0,0, 10), 00026 _corr01("Corr01", 0.0, -1.0, 1.0), 00027 _corr02("Corr02", 0.0, -1.0, 1.0), 00028 _corr12("Corr12", 0.0, -1.0, 1.0) 00029 {} 00030 00031 TrivariateGaussian::~TrivariateGaussian() { 00032 } 00033 00034 TrivariateGaussian::TrivariateGaussian(const TrivariateGaussian & right): 00035 _mean0(right._mean0), 00036 _mean1(right._mean1), 00037 _mean2(right._mean2), 00038 _sigma0(right._sigma0), 00039 _sigma1(right._sigma1), 00040 _sigma2(right._sigma2), 00041 _corr01(right._corr01), 00042 _corr02(right._corr02), 00043 _corr12(right._corr12) 00044 { 00045 } 00046 00047 00048 double TrivariateGaussian::operator() (const Argument & a) const { 00049 assert (a.dimension()==3); 00050 double x = a[0]; 00051 double y = a[1]; 00052 double z = a[2]; 00053 00054 00055 double x0 = _mean0.getValue(); 00056 double y0 = _mean1.getValue(); 00057 double z0 = _mean2.getValue(); 00058 00059 double dx = x-x0; 00060 double dy = y-y0; 00061 double dz = z-z0; 00062 00063 double sx = _sigma0.getValue(); 00064 double sy = _sigma1.getValue(); 00065 double sz = _sigma2.getValue(); 00066 00067 00068 double sxs = sx*sx; 00069 double sys = sy*sy; 00070 double szs = sz*sz; 00071 00072 00073 double rho1 = _corr01.getValue(); 00074 double rho2 = _corr12.getValue(); 00075 double rho3 = _corr02.getValue(); 00076 00077 double dt = (1.0+rho1*rho2*rho3-rho1*rho1-rho2*rho2-rho3*rho3); 00078 double tmp1 ,tmp2; 00079 00080 tmp1= 1.0/((2*M_PI)*sqrt(2*M_PI)*sx*sy*sz*sqrt(dt)); 00081 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)); 00082 00083 00084 return tmp1*tmp2; 00085 00086 00087 } 00088 00089 Parameter & TrivariateGaussian::mean0() { 00090 return _mean0; 00091 } 00092 00093 Parameter & TrivariateGaussian::sigma0() { 00094 return _sigma0; 00095 } 00096 00097 const Parameter & TrivariateGaussian::mean0() const { 00098 return _mean0; 00099 } 00100 00101 const Parameter & TrivariateGaussian::sigma0() const { 00102 return _sigma0; 00103 } 00104 00105 Parameter & TrivariateGaussian::mean1() { 00106 return _mean1; 00107 } 00108 00109 Parameter & TrivariateGaussian::sigma1() { 00110 return _sigma1; 00111 } 00112 00113 const Parameter & TrivariateGaussian::mean1() const { 00114 return _mean1; 00115 } 00116 00117 const Parameter & TrivariateGaussian::sigma1() const { 00118 return _sigma1; 00119 } 00120 00121 Parameter & TrivariateGaussian::mean2() { 00122 return _mean2; 00123 } 00124 00125 00126 Parameter & TrivariateGaussian::sigma2() { 00127 return _sigma2; 00128 } 00129 00130 const Parameter & TrivariateGaussian::mean2() const { 00131 return _mean2; 00132 } 00133 00134 const Parameter & TrivariateGaussian::sigma2() const { 00135 return _sigma2; 00136 } 00137 00138 00139 00140 Parameter & TrivariateGaussian::corr01() { 00141 return _corr01; 00142 } 00143 00144 const Parameter & TrivariateGaussian::corr01() const { 00145 return _corr01; 00146 } 00147 00148 Parameter & TrivariateGaussian::corr02() { 00149 return _corr02; 00150 } 00151 00152 const Parameter & TrivariateGaussian::corr02() const { 00153 return _corr02; 00154 } 00155 00156 Parameter & TrivariateGaussian::corr12() { 00157 return _corr12; 00158 } 00159 00160 const Parameter & TrivariateGaussian::corr12() const { 00161 return _corr12; 00162 } 00163 00164 00165 unsigned int TrivariateGaussian::dimensionality() const { 00166 return 3; 00167 } 00168 00169 double TrivariateGaussian::operator ()(double x) const 00170 { 00171 std::cerr 00172 << "Warning. trivariate Gaussian called with scalar argument" 00173 << std::endl; 00174 assert(0); 00175 return 0; 00176 } 00177 00178 } // namespace Genfun