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