CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 #include "CLHEP/GenericFunctions/VoigtProfile.hh" 00002 #include "CLHEP/GenericFunctions/Variable.hh" 00003 #include <cmath> 00004 #include <complex> 00005 00006 #if (defined __STRICT_ANSI__) || (defined _WIN32) 00007 #ifndef M_PI 00008 #define M_PI 3.14159265358979323846 00009 #endif // M_PI 00010 #endif // __STRICT_ANSI__ 00011 00012 using namespace std; 00013 00014 inline double Pow(double x,int n) { 00015 double val=1; 00016 for(int i=0;i<n;i++){ 00017 val=val*x; 00018 } 00019 return val; 00020 } 00021 00022 inline std::complex<double> nwwerf(std::complex<double> z) { 00023 std::complex<double> zh,r[38],s,t,v; 00024 00025 const double z1 = 1; 00026 const double hf = z1/2; 00027 const double z10 = 10; 00028 const double c1 = 74/z10; 00029 const double c2 = 83/z10; 00030 const double c3 = z10/32; 00031 const double c4 = 16/z10; 00032 const double c = 1.12837916709551257; 00033 const double p = Pow(2.0*c4,33); 00034 00035 double x=z.real(); 00036 double y=z.imag(); 00037 double xa=(x >= 0) ? x : -x; 00038 double ya=(y >= 0) ? y : -y; 00039 if(ya < c1 && xa < c2){ 00040 zh = std::complex<double>(ya+c4,xa); 00041 r[37]= std::complex<double>(0,0); 00042 // do 1 n = 36,1,-1 00043 for(int n = 36; n>0; n--){ 00044 t=zh+double(n)*std::conj(r[n+1]); 00045 r[n]=hf*t/std::norm(t); 00046 } 00047 double xl=p; 00048 s=std::complex<double>(0,0); 00049 // do 2 n = 33,1,-1 00050 for(int k=33; k>0; k--){ 00051 xl=c3*xl; 00052 s=r[k]*(s+xl); 00053 } 00054 v=c*s; 00055 } 00056 else{ 00057 zh=std::complex<double>(ya,xa); 00058 r[1]=std::complex<double>(0,0); 00059 // do 3 n = 9,1,-1 00060 for(int n=9;n>0;n--){ 00061 t=zh+double(n)*std::conj(r[1]); 00062 r[1]=hf*t/std::norm(t); 00063 } 00064 v=c*r[1]; 00065 } 00066 if(ya == 0) v= std::complex<double>(exp(-xa*xa),v.imag()); 00067 if(y < 0) { 00068 v=2.0*std::exp(std::complex<double>(-xa,-ya)*std::complex<double>(xa,ya))-v; 00069 if(x > 0) v=std::conj(v); 00070 } 00071 else{ 00072 if(x < 0) v=std::conj(v); 00073 } 00074 return v; 00075 } 00076 00077 00078 00079 namespace Genfun { 00080 FUNCTION_OBJECT_IMP(VoigtProfile) 00081 00082 00083 VoigtProfile::VoigtProfile(): 00084 _mass("mass", 50, 10, 90), 00085 _width ("width", 5, 0, 100), 00086 _sigma("sigma", 5, 0, 100) 00087 {} 00088 00089 VoigtProfile::VoigtProfile(const VoigtProfile & right): 00090 AbsFunction(), 00091 _mass(right._mass), 00092 _width (right._width), 00093 _sigma(right._sigma) 00094 { 00095 } 00096 00097 VoigtProfile::~VoigtProfile() { 00098 } 00099 00100 double VoigtProfile::operator() (double x) const { 00101 double M=_mass.getValue(); 00102 double G=_width.getValue()/2.0; 00103 double s=_sigma.getValue(); 00104 static const double sqrt2=sqrt(2.0); 00105 static const double sqrt2PI=sqrt(2.0*M_PI); 00106 static const std::complex<double> I(0,1); 00107 std::complex<double> z = ((x-M) + I*G)/sqrt2/s; 00108 00109 double f=nwwerf(z).real()/s/sqrt2PI; 00110 00111 return f; 00112 00113 } 00114 00115 Parameter & VoigtProfile::mass() { 00116 return _mass; 00117 } 00118 00119 00120 Parameter & VoigtProfile::width() { 00121 return _width; 00122 } 00123 00124 Parameter & VoigtProfile::sigma() { 00125 return _sigma; 00126 } 00127 00128 00129 } // namespace Genfun