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