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

AnalyticConvolution.cc
Go to the documentation of this file.
1// -*- C++ -*-
2// $Id: AnalyticConvolution.cc,v 1.8 2010/07/22 21:55:10 garren Exp $
6#include <cmath> // for isfinite
7#if (defined _WIN32)
8#include <float.h> // Visual C++ _finite
9#endif
10namespace Genfun {
11FUNCTION_OBJECT_IMP(AnalyticConvolution)
12
14 _lifetime ("Lifetime", 1.0, 0.0), // Bounded from below by zero, by default
15 _frequency("Frequency", 0.0, 0.0), // Bounded from below by zero, by default
16 _sigma ("Sigma", 1.0, 0.0), // Bounded from below by zero, by default
17 _offset ("Offset", 0.0), // Offset is unbounded
18 _type(type)
19{
20}
21
23 AbsFunction(right),
24 _lifetime (right._lifetime),
25 _frequency (right._frequency),
26 _sigma (right._sigma),
27 _offset (right._offset),
28 _type(right._type)
29{
30}
31
33}
34
36 return _frequency;
37}
38
40 return _frequency;
41}
42
44 return _lifetime;
45}
46
48 return _lifetime;
49}
50
52 return _sigma;
53}
54
56 return _sigma;
57}
58
60 return _offset;
61}
62
64 return _offset;
65}
66double AnalyticConvolution::operator() (double argument) const {
67 // Fetch the paramters. This operator does not convolve numerically.
68 static const double sqrtTwo = sqrt(2.0);
69 double xsigma = _sigma.getValue();
70 double tau = _lifetime.getValue();
71 double xoffset = _offset.getValue();
72 double x = argument-xoffset;
73 double freq = _frequency.getValue();
74
75 // smeared exponential an its asymmetry.
76 double expG=0.0, asymm=0.0;
77
78 if (_type==SMEARED_NEG_EXP) {
79 expG = exp((xsigma*xsigma +2*tau*(/*xoffset*/x))/(2.0*tau*tau)) *
80 erfc((xsigma*xsigma+tau*(/*xoffset*/x))/(sqrtTwo*xsigma*tau))/(2.0*tau);
81#if (defined _WIN32)
82 if (!_finite(expG)) {
83 expG=0.0;
84 }
85#else
86 if (!std::isfinite(expG)) {
87 expG=0.0;
88 }
89#endif
90 return expG;
91 }
92 else {
93 expG = exp((xsigma*xsigma +2*tau*(/*xoffset*/-x))/(2.0*tau*tau)) *
94 erfc((xsigma*xsigma+tau*(/*xoffset*/-x))/(sqrtTwo*xsigma*tau))/(2.0*tau);
95 }
96
97 // Both sign distribution=> return smeared exponential:
98 if (_type==SMEARED_EXP) {
99#if (defined _WIN32)
100 if (!_finite(expG)) {
101 expG=0.0;
102 }
103#else
104 if (!std::isfinite(expG)) {
105 expG=0.0;
106 }
107#endif
108 return expG;
109 }
110
111
112 // Calcualtion of asymmetry:
113
114 // First, if the mixing frequency is too high compared with the lifetime, we
115 // cannot see this oscillation. We abandon the complicated approach and just do
116 // this instead:
117 if (xsigma>6.0*tau) {
118 asymm = expG*(1/(1+tau*tau*freq*freq));
119 }
120 else if (xsigma==0.0) {
121 if (_type==SMEARED_COS_EXP|| _type==MIXED || _type==UNMIXED ) {
122 if (x>=0) asymm= (expG*cos(freq*x));
123 }
124 else if (_type==SMEARED_SIN_EXP) {
125 if (x>=0) asymm= (expG*sin(freq*x));
126 }
127 }
128 else {
129 std::complex<double> z(freq*xsigma/sqrtTwo, (xsigma/tau-x/xsigma)/sqrtTwo);
130 if (x<0) {
131 if (_type==SMEARED_COS_EXP|| _type==MIXED || _type==UNMIXED ) {
132 asymm= 2.0*nwwerf(z).real()/tau/4.0*exp(-x*x/2.0/xsigma/xsigma);
133 }
134 else if (_type==SMEARED_SIN_EXP) {
135 asymm= 2.0*nwwerf(z).imag()/tau/4.0*exp(-x*x/2.0/xsigma/xsigma);
136 }
137 }
138 else {
139 if (_type==SMEARED_COS_EXP||_type==MIXED || _type==UNMIXED) {
140 asymm= -2.0*nwwerf(std::conj(z)).real()/tau/4*exp(-x*x/2.0/xsigma/xsigma) +
141 exp(xsigma*xsigma/2 *(1/tau/tau - freq*freq) - x/tau)*(1./tau)*cos(freq*x - freq/tau*xsigma*xsigma);
142 }
143 else if (_type==SMEARED_SIN_EXP) {
144 asymm= +2.0*nwwerf(std::conj(z)).imag()/tau/4*exp(-x*x/2.0/xsigma/xsigma) +
145 exp(xsigma*xsigma/2 *(1/tau/tau - freq*freq) - x/tau)*(1./tau)*sin(freq*x - freq/tau*xsigma*xsigma);
146 }
147 }
148 }
149
150 // Return either MIXED, UNMIXED, or ASYMMETRY function.
151 if (_type==UNMIXED) {
152 double retVal = (expG+asymm)/2.0;
153 if (retVal<0)
154 std::cerr
155 << "Warning in AnalyticConvolution: negative probablity"
156 << std::endl;
157 if (retVal<0)
158 std::cerr
159 << xsigma << ' ' << tau << ' ' << xoffset << ' '
160 << freq << ' '<< argument
161 << std::endl;
162 if (retVal<0)
163 std::cerr << retVal << std::endl;
164 return retVal;
165 }
166 else if (_type==MIXED) {
167 double retVal = (expG-asymm)/2.0;
168 if (retVal<0)
169 std::cerr
170 << "Warning in AnalyticConvolution: negative probablity"
171 << std::endl;
172 if (retVal<0)
173 std::cerr
174 << xsigma << ' ' << tau << ' ' << xoffset << ' '
175 << freq << ' ' << argument
176 << std::endl;
177 if (retVal<0)
178 std::cerr << retVal << std::endl;
179 return retVal;
180 }
181 else if (_type==SMEARED_COS_EXP || _type==SMEARED_SIN_EXP) {
182 return asymm;
183 }
184 else {
185 std::cerr
186 << "Unknown sign parity. State is not allowed"
187 << std::endl;
188 exit(0);
189 return 0.0;
190 }
191
192}
193
194double AnalyticConvolution::erfc(double x) const {
195// This is accurate to 7 places.
196// See Numerical Recipes P 221
197 double t, z, ans;
198 z = (x < 0) ? -x : x;
199 t = 1.0/(1.0+.5*z);
200 ans = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
201 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
202 t*(-0.82215223+t*0.17087277 ))) ))) )));
203 if ( x < 0 ) ans = 2.0 - ans;
204 return ans;
205}
206
207double AnalyticConvolution::pow(double x,int n) const {
208 double val=1;
209 for(int i=0;i<n;i++){
210 val=val*x;
211 }
212 return val;
213}
214
215std::complex<double> AnalyticConvolution::nwwerf(std::complex<double> z) const {
216 std::complex<double> zh,r[38],s,t,v;
217
218 const double z1 = 1;
219 const double hf = z1/2;
220 const double z10 = 10;
221 const double c1 = 74/z10;
222 const double c2 = 83/z10;
223 const double c3 = z10/32;
224 const double c4 = 16/z10;
225 const double c = 1.12837916709551257;
226 const double p = pow(2.0*c4,33);
227
228 double x=z.real();
229 double y=z.imag();
230 double xa=(x >= 0) ? x : -x;
231 double ya=(y >= 0) ? y : -y;
232 if(ya < c1 && xa < c2){
233 zh = std::complex<double>(ya+c4,xa);
234 r[37]= std::complex<double>(0,0);
235 // do 1 n = 36,1,-1
236 for(int n = 36; n>0; n--){
237 t=zh+double(n)*std::conj(r[n+1]);
238 r[n]=hf*t/std::norm(t);
239 }
240 double xl=p;
241 s=std::complex<double>(0,0);
242 // do 2 n = 33,1,-1
243 for(int k=33; k>0; k--){
244 xl=c3*xl;
245 s=r[k]*(s+xl);
246 }
247 v=c*s;
248 }
249 else{
250 zh=std::complex<double>(ya,xa);
251 r[1]=std::complex<double>(0,0);
252 // do 3 n = 9,1,-1
253 for(int n=9;n>0;n--){
254 t=zh+double(n)*std::conj(r[1]);
255 r[1]=hf*t/std::norm(t);
256 }
257 v=c*r[1];
258 }
259 if(ya == 0) v= std::complex<double>(exp(-xa*xa),v.imag());
260 if(y < 0) {
261 v=2.0*std::exp(std::complex<double>(-xa,-ya)*std::complex<double>(xa,ya))-v;
262 if(x > 0) v=std::conj(v);
263 }
264 else{
265 if(x < 0) v=std::conj(v);
266 }
267 return v;
268}
269} // namespace Genfun
#define FUNCTION_OBJECT_IMP(classname)
virtual double operator()(double argument) const
AnalyticConvolution(Type=SMEARED_EXP)
virtual double getValue() const
Definition: Parameter.cc:27
#define double(obj)
Definition: excDblThrow.cc:32
#define exit(x)
double norm(const HepGenMatrix &m)
Definition: GenMatrix.cc:57