CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: PuncturedSmearedExp.cc,v 1.2.4.1 2004/04/29 00:20:36 garren Exp $ 00003 #include "CLHEP/GenericFunctions/PuncturedSmearedExp.hh" 00004 #include <sstream> 00005 #include <cmath> 00006 namespace Genfun { 00007 FUNCTION_OBJECT_IMP(PuncturedSmearedExp) 00008 00009 PuncturedSmearedExp::PuncturedSmearedExp() : 00010 _lifetime ("Lifetime", 1.0, 0.0), // Bounded from below by zero, by default 00011 _sigma ("Sigma", 1.0, 0.0) // Bounded from below by zero, by default 00012 { 00013 } 00014 00015 PuncturedSmearedExp::PuncturedSmearedExp(const PuncturedSmearedExp & right): 00016 _lifetime (right._lifetime), 00017 _sigma (right._sigma), 00018 _punctures (right._punctures) 00019 { 00020 } 00021 00022 void PuncturedSmearedExp::puncture(double min, double max) { 00023 std::ostringstream mn, mx; 00024 mn << "Min_" << _punctures.size()/2; 00025 mx << "Max_" << _punctures.size()/2; 00026 _punctures.push_back(Parameter(mn.str(), min, 0.0, 10.0)); 00027 _punctures.push_back(Parameter(mx.str(), max, 0.0, 10.0)); 00028 } 00029 00030 00031 Parameter & PuncturedSmearedExp::min(unsigned int i) { 00032 return _punctures[2*i]; 00033 } 00034 00035 const Parameter & PuncturedSmearedExp::min(unsigned int i) const { 00036 return _punctures[2*i]; 00037 } 00038 00039 00040 Parameter & PuncturedSmearedExp::max(unsigned int i) { 00041 return _punctures[2*i+1]; 00042 00043 } 00044 00045 const Parameter & PuncturedSmearedExp::max(unsigned int i) const { 00046 return _punctures[2*i+1]; 00047 } 00048 00049 00050 PuncturedSmearedExp::~PuncturedSmearedExp() { 00051 } 00052 00053 Parameter & PuncturedSmearedExp::lifetime() { 00054 return _lifetime; 00055 } 00056 00057 const Parameter & PuncturedSmearedExp::lifetime() const { 00058 return _lifetime; 00059 } 00060 00061 Parameter & PuncturedSmearedExp::sigma() { 00062 return _sigma; 00063 } 00064 00065 const Parameter & PuncturedSmearedExp::sigma() const { 00066 return _sigma; 00067 } 00068 00069 00070 double PuncturedSmearedExp::operator() (double argument) const { 00071 // Fetch the paramters. This operator does not convolve numerically. 00072 static const double sqrtTwo = sqrt(2.0); 00073 00074 double sigma = _sigma.getValue(); 00075 double tau = _lifetime.getValue(); 00076 double x = argument; 00077 00078 // Copy: 00079 std::vector<double> punctures(_punctures.size()); 00080 for (size_t i=0;i<_punctures.size();i++) punctures[i]=_punctures[i].getValue(); 00081 00082 // Overlap elimination: 00083 bool overlap=true; 00084 00085 while (overlap) { 00086 00087 overlap=false; 00088 00089 for (size_t i=0;i<punctures.size()/2;i++) { 00090 std::sort(punctures.begin()+2*i, punctures.begin()+2*i+2); 00091 double min1=punctures[2*i]; 00092 double max1=punctures[2*i+1]; 00093 for (size_t j=i+1;j<punctures.size()/2;j++) { 00094 std::sort(punctures.begin()+2*j, punctures.begin()+2*j+2); 00095 double min2=punctures[2*j]; 00096 double max2=punctures[2*j+1]; 00097 if ((min2>min1 && max1>min2) || (min1>min2 && max2<min1)) { 00098 punctures[2*i] =std::min(min1,min2); 00099 punctures[2*i+1]=std::max(max1,max2); 00100 std::vector<double>::iterator t0=punctures.begin()+2*j, t1=t0+2; 00101 punctures.erase(t0,t1); 00102 overlap=true; 00103 break; 00104 } 00105 } 00106 if (overlap) break; 00107 } 00108 } 00109 00110 double expG=0,norm=0; 00111 for (size_t i=0;i<punctures.size()/2;i++) { 00112 00113 double a = punctures[2*i]; 00114 double b = punctures[2*i+1]; 00115 00116 double alpha = (a/sigma + sigma/tau)/sqrtTwo; 00117 double beta = (b/sigma + sigma/tau)/sqrtTwo; 00118 double delta = 1/sqrtTwo/sigma; 00119 00120 00121 norm += 2*tau*exp(1/(4*delta*delta*tau*tau))*(exp(-alpha/(delta*tau)) - exp(-beta/(delta*tau))); 00122 00123 expG += ((erfc(alpha-delta*x) - erfc(beta-delta*x))*exp(-x/tau) ); 00124 00125 00126 } 00127 00128 // protection: 00129 if (norm==0) return norm; 00130 00131 return expG/norm; 00132 } 00133 00134 double PuncturedSmearedExp::erfc(double x) const { 00135 // This is accurate to 7 places. 00136 // See Numerical Recipes P 221 00137 double t, z, ans; 00138 z = (x < 0) ? -x : x; 00139 t = 1.0/(1.0+.5*z); 00140 ans = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+ 00141 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+ 00142 t*(-0.82215223+t*0.17087277 ))) ))) ))); 00143 if ( x < 0 ) ans = 2.0 - ans; 00144 return ans; 00145 } 00146 00147 double PuncturedSmearedExp::pow(double x,int n) const { 00148 double val=1; 00149 for(int i=0;i<n;i++){ 00150 val=val*x; 00151 } 00152 return val; 00153 } 00154 00155 } // namespace Genfun 00156 00157