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

VoigtProfile.cc
Go to the documentation of this file.
3 #include <cmath>
4 #include <complex>
5 
6 #if (defined __STRICT_ANSI__) || (defined _WIN32)
7 #ifndef M_PI
8 #define M_PI 3.14159265358979323846
9 #endif // M_PI
10 #endif // __STRICT_ANSI__
11 
12 using namespace std;
13 
14 inline double Pow(double x,int n) {
15  double val=1;
16  for(int i=0;i<n;i++){
17  val=val*x;
18  }
19  return val;
20 }
21 
22 inline std::complex<double> nwwerf(std::complex<double> z) {
23  std::complex<double> zh,r[38],s,t,v;
24 
25  const double z1 = 1;
26  const double hf = z1/2;
27  const double z10 = 10;
28  const double c1 = 74/z10;
29  const double c2 = 83/z10;
30  const double c3 = z10/32;
31  const double c4 = 16/z10;
32  const double c = 1.12837916709551257;
33  const double p = Pow(2.0*c4,33);
34 
35  double x=z.real();
36  double y=z.imag();
37  double xa=(x >= 0) ? x : -x;
38  double ya=(y >= 0) ? y : -y;
39  if(ya < c1 && xa < c2){
40  zh = std::complex<double>(ya+c4,xa);
41  r[37]= std::complex<double>(0,0);
42  // do 1 n = 36,1,-1
43  for(int n = 36; n>0; n--){
44  t=zh+double(n)*std::conj(r[n+1]);
45  r[n]=hf*t/std::norm(t);
46  }
47  double xl=p;
48  s=std::complex<double>(0,0);
49  // do 2 n = 33,1,-1
50  for(int k=33; k>0; k--){
51  xl=c3*xl;
52  s=r[k]*(s+xl);
53  }
54  v=c*s;
55  }
56  else{
57  zh=std::complex<double>(ya,xa);
58  r[1]=std::complex<double>(0,0);
59  // do 3 n = 9,1,-1
60  for(int n=9;n>0;n--){
61  t=zh+double(n)*std::conj(r[1]);
62  r[1]=hf*t/std::norm(t);
63  }
64  v=c*r[1];
65  }
66  if(ya == 0) v= std::complex<double>(exp(-xa*xa),v.imag());
67  if(y < 0) {
68  v=2.0*std::exp(std::complex<double>(-xa,-ya)*std::complex<double>(xa,ya))-v;
69  if(x > 0) v=std::conj(v);
70  }
71  else{
72  if(x < 0) v=std::conj(v);
73  }
74  return v;
75 }
76 
77 
78 
79 namespace Genfun {
80 FUNCTION_OBJECT_IMP(VoigtProfile)
81 
82 
84  _mass("mass", 50, 10, 90),
85  _width ("width", 5, 0, 100),
86  _sigma("sigma", 5, 0, 100)
87 {}
88 
89  VoigtProfile::VoigtProfile(const VoigtProfile & right):
90  AbsFunction(),
91  _mass(right._mass),
92  _width (right._width),
93  _sigma(right._sigma)
94 {
95 }
96 
98 }
99 
100 double VoigtProfile::operator() (double x) const {
101  double M=_mass.getValue();
102  double G=_width.getValue()/2.0;
103  double s=_sigma.getValue();
104  static const double sqrt2=sqrt(2.0);
105  static const double sqrt2PI=sqrt(2.0*M_PI);
106  static const std::complex<double> I(0,1);
107  std::complex<double> z = ((x-M) + I*G)/sqrt2/s;
108 
109  double f=nwwerf(z).real()/s/sqrt2PI;
110 
111  return f;
112 
113 }
114 
116  return _mass;
117 }
118 
119 
121  return _width;
122 }
123 
125  return _sigma;
126 }
127 
128 
129 } // namespace Genfun