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