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