CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: AnalyticConvolution.cc,v 1.8 2010/07/22 21:55:10 garren Exp $ 00003 #include "CLHEP/GenericFunctions/AnalyticConvolution.hh" 00004 #include "CLHEP/GenericFunctions/Gaussian.hh" 00005 #include "CLHEP/GenericFunctions/Exponential.hh" 00006 #include <cmath> // for isfinite 00007 #if (defined _WIN32) 00008 #include <float.h> // Visual C++ _finite 00009 #endif 00010 namespace Genfun { 00011 FUNCTION_OBJECT_IMP(AnalyticConvolution) 00012 00013 AnalyticConvolution::AnalyticConvolution(AnalyticConvolution::Type type) : 00014 _lifetime ("Lifetime", 1.0, 0.0), // Bounded from below by zero, by default 00015 _frequency("Frequency", 0.0, 0.0), // Bounded from below by zero, by default 00016 _sigma ("Sigma", 1.0, 0.0), // Bounded from below by zero, by default 00017 _offset ("Offset", 0.0), // Offset is unbounded 00018 _type(type) 00019 { 00020 } 00021 00022 AnalyticConvolution::AnalyticConvolution(const AnalyticConvolution & right): 00023 AbsFunction(right), 00024 _lifetime (right._lifetime), 00025 _frequency (right._frequency), 00026 _sigma (right._sigma), 00027 _offset (right._offset), 00028 _type(right._type) 00029 { 00030 } 00031 00032 AnalyticConvolution::~AnalyticConvolution() { 00033 } 00034 00035 Parameter & AnalyticConvolution::frequency() { 00036 return _frequency; 00037 } 00038 00039 const Parameter & AnalyticConvolution::frequency() const { 00040 return _frequency; 00041 } 00042 00043 Parameter & AnalyticConvolution::lifetime() { 00044 return _lifetime; 00045 } 00046 00047 const Parameter & AnalyticConvolution::lifetime() const { 00048 return _lifetime; 00049 } 00050 00051 Parameter & AnalyticConvolution::sigma() { 00052 return _sigma; 00053 } 00054 00055 const Parameter & AnalyticConvolution::sigma() const { 00056 return _sigma; 00057 } 00058 00059 Parameter & AnalyticConvolution::offset() { 00060 return _offset; 00061 } 00062 00063 const Parameter & AnalyticConvolution::offset() const { 00064 return _offset; 00065 } 00066 double AnalyticConvolution::operator() (double argument) const { 00067 // Fetch the paramters. This operator does not convolve numerically. 00068 static const double sqrtTwo = sqrt(2.0); 00069 double xsigma = _sigma.getValue(); 00070 double tau = _lifetime.getValue(); 00071 double xoffset = _offset.getValue(); 00072 double x = argument-xoffset; 00073 double freq = _frequency.getValue(); 00074 00075 // smeared exponential an its asymmetry. 00076 double expG=0.0, asymm=0.0; 00077 00078 if (_type==SMEARED_NEG_EXP) { 00079 expG = exp((xsigma*xsigma +2*tau*(/*xoffset*/x))/(2.0*tau*tau)) * 00080 erfc((xsigma*xsigma+tau*(/*xoffset*/x))/(sqrtTwo*xsigma*tau))/(2.0*tau); 00081 #if (defined _WIN32) 00082 if (!_finite(expG)) { 00083 expG=0.0; 00084 } 00085 #else 00086 if (!std::isfinite(expG)) { 00087 expG=0.0; 00088 } 00089 #endif 00090 return expG; 00091 } 00092 else { 00093 expG = exp((xsigma*xsigma +2*tau*(/*xoffset*/-x))/(2.0*tau*tau)) * 00094 erfc((xsigma*xsigma+tau*(/*xoffset*/-x))/(sqrtTwo*xsigma*tau))/(2.0*tau); 00095 } 00096 00097 // Both sign distribution=> return smeared exponential: 00098 if (_type==SMEARED_EXP) { 00099 #if (defined _WIN32) 00100 if (!_finite(expG)) { 00101 expG=0.0; 00102 } 00103 #else 00104 if (!std::isfinite(expG)) { 00105 expG=0.0; 00106 } 00107 #endif 00108 return expG; 00109 } 00110 00111 00112 // Calcualtion of asymmetry: 00113 00114 // First, if the mixing frequency is too high compared with the lifetime, we 00115 // cannot see this oscillation. We abandon the complicated approach and just do 00116 // this instead: 00117 if (xsigma>6.0*tau) { 00118 asymm = expG*(1/(1+tau*tau*freq*freq)); 00119 } 00120 else if (xsigma==0.0) { 00121 if (_type==SMEARED_COS_EXP|| _type==MIXED || _type==UNMIXED ) { 00122 if (x>=0) asymm= (expG*cos(freq*x)); 00123 } 00124 else if (_type==SMEARED_SIN_EXP) { 00125 if (x>=0) asymm= (expG*sin(freq*x)); 00126 } 00127 } 00128 else { 00129 std::complex<double> z(freq*xsigma/sqrtTwo, (xsigma/tau-x/xsigma)/sqrtTwo); 00130 if (x<0) { 00131 if (_type==SMEARED_COS_EXP|| _type==MIXED || _type==UNMIXED ) { 00132 asymm= 2.0*nwwerf(z).real()/tau/4.0*exp(-x*x/2.0/xsigma/xsigma); 00133 } 00134 else if (_type==SMEARED_SIN_EXP) { 00135 asymm= 2.0*nwwerf(z).imag()/tau/4.0*exp(-x*x/2.0/xsigma/xsigma); 00136 } 00137 } 00138 else { 00139 if (_type==SMEARED_COS_EXP||_type==MIXED || _type==UNMIXED) { 00140 asymm= -2.0*nwwerf(std::conj(z)).real()/tau/4*exp(-x*x/2.0/xsigma/xsigma) + 00141 exp(xsigma*xsigma/2 *(1/tau/tau - freq*freq) - x/tau)*(1./tau)*cos(freq*x - freq/tau*xsigma*xsigma); 00142 } 00143 else if (_type==SMEARED_SIN_EXP) { 00144 asymm= +2.0*nwwerf(std::conj(z)).imag()/tau/4*exp(-x*x/2.0/xsigma/xsigma) + 00145 exp(xsigma*xsigma/2 *(1/tau/tau - freq*freq) - x/tau)*(1./tau)*sin(freq*x - freq/tau*xsigma*xsigma); 00146 } 00147 } 00148 } 00149 00150 // Return either MIXED, UNMIXED, or ASYMMETRY function. 00151 if (_type==UNMIXED) { 00152 double retVal = (expG+asymm)/2.0; 00153 if (retVal<0) 00154 std::cerr 00155 << "Warning in AnalyticConvolution: negative probablity" 00156 << std::endl; 00157 if (retVal<0) 00158 std::cerr 00159 << xsigma << ' ' << tau << ' ' << xoffset << ' ' 00160 << freq << ' '<< argument 00161 << std::endl; 00162 if (retVal<0) 00163 std::cerr << retVal << std::endl; 00164 return retVal; 00165 } 00166 else if (_type==MIXED) { 00167 double retVal = (expG-asymm)/2.0; 00168 if (retVal<0) 00169 std::cerr 00170 << "Warning in AnalyticConvolution: negative probablity" 00171 << std::endl; 00172 if (retVal<0) 00173 std::cerr 00174 << xsigma << ' ' << tau << ' ' << xoffset << ' ' 00175 << freq << ' ' << argument 00176 << std::endl; 00177 if (retVal<0) 00178 std::cerr << retVal << std::endl; 00179 return retVal; 00180 } 00181 else if (_type==SMEARED_COS_EXP || _type==SMEARED_SIN_EXP) { 00182 return asymm; 00183 } 00184 else { 00185 std::cerr 00186 << "Unknown sign parity. State is not allowed" 00187 << std::endl; 00188 exit(0); 00189 return 0.0; 00190 } 00191 00192 } 00193 00194 double AnalyticConvolution::erfc(double x) const { 00195 // This is accurate to 7 places. 00196 // See Numerical Recipes P 221 00197 double t, z, ans; 00198 z = (x < 0) ? -x : x; 00199 t = 1.0/(1.0+.5*z); 00200 ans = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+ 00201 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+ 00202 t*(-0.82215223+t*0.17087277 ))) ))) ))); 00203 if ( x < 0 ) ans = 2.0 - ans; 00204 return ans; 00205 } 00206 00207 double AnalyticConvolution::pow(double x,int n) const { 00208 double val=1; 00209 for(int i=0;i<n;i++){ 00210 val=val*x; 00211 } 00212 return val; 00213 } 00214 00215 std::complex<double> AnalyticConvolution::nwwerf(std::complex<double> z) const { 00216 std::complex<double> zh,r[38],s,t,v; 00217 00218 const double z1 = 1; 00219 const double hf = z1/2; 00220 const double z10 = 10; 00221 const double c1 = 74/z10; 00222 const double c2 = 83/z10; 00223 const double c3 = z10/32; 00224 const double c4 = 16/z10; 00225 const double c = 1.12837916709551257; 00226 const double p = pow(2.0*c4,33); 00227 00228 double x=z.real(); 00229 double y=z.imag(); 00230 double xa=(x >= 0) ? x : -x; 00231 double ya=(y >= 0) ? y : -y; 00232 if(ya < c1 && xa < c2){ 00233 zh = std::complex<double>(ya+c4,xa); 00234 r[37]= std::complex<double>(0,0); 00235 // do 1 n = 36,1,-1 00236 for(int n = 36; n>0; n--){ 00237 t=zh+double(n)*std::conj(r[n+1]); 00238 r[n]=hf*t/std::norm(t); 00239 } 00240 double xl=p; 00241 s=std::complex<double>(0,0); 00242 // do 2 n = 33,1,-1 00243 for(int k=33; k>0; k--){ 00244 xl=c3*xl; 00245 s=r[k]*(s+xl); 00246 } 00247 v=c*s; 00248 } 00249 else{ 00250 zh=std::complex<double>(ya,xa); 00251 r[1]=std::complex<double>(0,0); 00252 // do 3 n = 9,1,-1 00253 for(int n=9;n>0;n--){ 00254 t=zh+double(n)*std::conj(r[1]); 00255 r[1]=hf*t/std::norm(t); 00256 } 00257 v=c*r[1]; 00258 } 00259 if(ya == 0) v= std::complex<double>(exp(-xa*xa),v.imag()); 00260 if(y < 0) { 00261 v=2.0*std::exp(std::complex<double>(-xa,-ya)*std::complex<double>(xa,ya))-v; 00262 if(x > 0) v=std::conj(v); 00263 } 00264 else{ 00265 if(x < 0) v=std::conj(v); 00266 } 00267 return v; 00268 } 00269 } // namespace Genfun