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  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 }
66 double 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 
194 double 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 
207 double 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 
215 std::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 double(obj)
Definition: excDblThrow.cc:32
double norm(const HepGenMatrix &m)
Definition: GenMatrix.cc:57
virtual double operator()(double argument) const
AnalyticConvolution(Type=SMEARED_EXP)
#define exit(x)
virtual double getValue() const
Definition: Parameter.cc:27
#define FUNCTION_OBJECT_IMP(classname)