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
Genfun::Parameter::getValue
virtual double getValue() const
Definition: Parameter.cc:27
double
#define double(obj)
Definition: excDblThrow.cc:32
Genfun::VoigtProfile
Definition: CLHEP/GenericFunctions/VoigtProfile.hh:20
Genfun::AbsFunction
Definition: CLHEP/GenericFunctions/AbsFunction.hh:48
Genfun::VoigtProfile::sigma
Parameter & sigma()
Definition: VoigtProfile.cc:124
Pow
double Pow(double x, int n)
Definition: VoigtProfile.cc:14
VoigtProfile.hh
nwwerf
std::complex< double > nwwerf(std::complex< double > z)
Definition: VoigtProfile.cc:22
CLHEP::detail::n
n
Definition: Ranlux64Engine.cc:85
Genfun::VoigtProfile::width
Parameter & width()
Definition: VoigtProfile.cc:120
f
void f(void g())
Definition: excDblThrow.cc:38
Variable.hh
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
Genfun::VoigtProfile::mass
Parameter & mass()
Definition: VoigtProfile.cc:115
s
Methods applicble to containers of as in std::list< LorentzVector > s
Definition: keyMergeIssues.doc:328
i
long i
Definition: JamesRandomSeeding.txt:27
Genfun::Parameter
Definition: CLHEP/GenericFunctions/Parameter.hh:35
x
any side effects of that construction would occur twice The semantics of throw x
Definition: whyZMthrowRethrows.txt:37
k
long k
Definition: JamesRandomSeeding.txt:29
FUNCTION_OBJECT_IMP
#define FUNCTION_OBJECT_IMP(classname)
Definition: CLHEP/GenericFunctions/AbsFunction.hh:156
Genfun::VoigtProfile::~VoigtProfile
virtual ~VoigtProfile()
Definition: VoigtProfile.cc:97
CLHEP::norm
double norm(const HepGenMatrix &m)
Definition: GenMatrix.cc:57
Genfun
Definition: CLHEP/GenericFunctions/Abs.hh:14
Genfun::VoigtProfile::operator()
virtual double operator()(double argument) const
Definition: VoigtProfile.cc:100