CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: TrivariateGaussian.cc,v 1.8 2010/06/16 18:22:01 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 AbsFunction(right), 00036 _mean0(right._mean0), 00037 _mean1(right._mean1), 00038 _mean2(right._mean2), 00039 _sigma0(right._sigma0), 00040 _sigma1(right._sigma1), 00041 _sigma2(right._sigma2), 00042 _corr01(right._corr01), 00043 _corr02(right._corr02), 00044 _corr12(right._corr12) 00045 { 00046 } 00047 00048 00049 double TrivariateGaussian::operator() (const Argument & a) const { 00050 assert (a.dimension()==3); 00051 double x = a[0]; 00052 double y = a[1]; 00053 double z = a[2]; 00054 00055 00056 double x0 = _mean0.getValue(); 00057 double y0 = _mean1.getValue(); 00058 double z0 = _mean2.getValue(); 00059 00060 double dx = x-x0; 00061 double dy = y-y0; 00062 double dz = z-z0; 00063 00064 double sx = _sigma0.getValue(); 00065 double sy = _sigma1.getValue(); 00066 double sz = _sigma2.getValue(); 00067 00068 00069 double sxs = sx*sx; 00070 double sys = sy*sy; 00071 double szs = sz*sz; 00072 00073 00074 double rho1 = _corr01.getValue(); 00075 double rho2 = _corr12.getValue(); 00076 double rho3 = _corr02.getValue(); 00077 00078 double dt = (1.0+rho1*rho2*rho3-rho1*rho1-rho2*rho2-rho3*rho3); 00079 double tmp1 ,tmp2; 00080 00081 tmp1= 1.0/((2*M_PI)*sqrt(2*M_PI)*sx*sy*sz*sqrt(dt)); 00082 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)); 00083 00084 00085 return tmp1*tmp2; 00086 00087 00088 } 00089 00090 Parameter & TrivariateGaussian::mean0() { 00091 return _mean0; 00092 } 00093 00094 Parameter & TrivariateGaussian::sigma0() { 00095 return _sigma0; 00096 } 00097 00098 const Parameter & TrivariateGaussian::mean0() const { 00099 return _mean0; 00100 } 00101 00102 const Parameter & TrivariateGaussian::sigma0() const { 00103 return _sigma0; 00104 } 00105 00106 Parameter & TrivariateGaussian::mean1() { 00107 return _mean1; 00108 } 00109 00110 Parameter & TrivariateGaussian::sigma1() { 00111 return _sigma1; 00112 } 00113 00114 const Parameter & TrivariateGaussian::mean1() const { 00115 return _mean1; 00116 } 00117 00118 const Parameter & TrivariateGaussian::sigma1() const { 00119 return _sigma1; 00120 } 00121 00122 Parameter & TrivariateGaussian::mean2() { 00123 return _mean2; 00124 } 00125 00126 00127 Parameter & TrivariateGaussian::sigma2() { 00128 return _sigma2; 00129 } 00130 00131 const Parameter & TrivariateGaussian::mean2() const { 00132 return _mean2; 00133 } 00134 00135 const Parameter & TrivariateGaussian::sigma2() const { 00136 return _sigma2; 00137 } 00138 00139 00140 00141 Parameter & TrivariateGaussian::corr01() { 00142 return _corr01; 00143 } 00144 00145 const Parameter & TrivariateGaussian::corr01() const { 00146 return _corr01; 00147 } 00148 00149 Parameter & TrivariateGaussian::corr02() { 00150 return _corr02; 00151 } 00152 00153 const Parameter & TrivariateGaussian::corr02() const { 00154 return _corr02; 00155 } 00156 00157 Parameter & TrivariateGaussian::corr12() { 00158 return _corr12; 00159 } 00160 00161 const Parameter & TrivariateGaussian::corr12() const { 00162 return _corr12; 00163 } 00164 00165 00166 unsigned int TrivariateGaussian::dimensionality() const { 00167 return 3; 00168 } 00169 00170 double TrivariateGaussian::operator ()(double) const 00171 { 00172 std::cerr 00173 << "Warning. trivariate Gaussian called with scalar argument" 00174 << std::endl; 00175 assert(0); 00176 return 0; 00177 } 00178 00179 } // namespace Genfun