CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

AnalyticConvolution.cc
Go to the documentation of this file.
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