14 _lifetime ("Lifetime", 1.0, 0.0),
15 _frequency("Frequency", 0.0, 0.0),
16 _sigma ("
Sigma", 1.0, 0.0),
17 _offset ("Offset", 0.0),
24 _lifetime (right._lifetime),
25 _frequency (right._frequency),
26 _sigma (right._sigma),
27 _offset (right._offset),
68 static const double sqrtTwo = sqrt(2.0);
72 double x = argument-xoffset;
76 double expG=0.0, asymm=0.0;
79 expG = exp((xsigma*xsigma +2*tau*(x))/(2.0*tau*tau)) *
80 erfc((xsigma*xsigma+tau*(x))/(sqrtTwo*xsigma*tau))/(2.0*tau);
86 if (!std::isfinite(expG)) {
93 expG = exp((xsigma*xsigma +2*tau*(-x))/(2.0*tau*tau)) *
94 erfc((xsigma*xsigma+tau*(-x))/(sqrtTwo*xsigma*tau))/(2.0*tau);
100 if (!_finite(expG)) {
104 if (!std::isfinite(expG)) {
117 if (xsigma>6.0*tau) {
118 asymm = expG*(1/(1+tau*tau*freq*freq));
120 else if (xsigma==0.0) {
122 if (x>=0) asymm= (expG*cos(freq*x));
125 if (x>=0) asymm= (expG*sin(freq*x));
129 std::complex<double> z(freq*xsigma/sqrtTwo, (xsigma/tau-x/xsigma)/sqrtTwo);
132 asymm= 2.0*nwwerf(z).real()/tau/4.0*exp(-x*x/2.0/xsigma/xsigma);
135 asymm= 2.0*nwwerf(z).imag()/tau/4.0*exp(-x*x/2.0/xsigma/xsigma);
140 asymm= -2.0*nwwerf(std::conj(z)).real()/tau/4*exp(-x*x/2.0/xsigma/xsigma) +
141 exp(xsigma*xsigma/2 *(1/tau/tau - freq*freq) - x/tau)*(1./tau)*cos(freq*x - freq/tau*xsigma*xsigma);
144 asymm= +2.0*nwwerf(std::conj(z)).imag()/tau/4*exp(-x*x/2.0/xsigma/xsigma) +
145 exp(xsigma*xsigma/2 *(1/tau/tau - freq*freq) - x/tau)*(1./tau)*sin(freq*x - freq/tau*xsigma*xsigma);
152 double retVal = (expG+asymm)/2.0;
155 <<
"Warning in AnalyticConvolution: negative probablity"
159 << xsigma <<
' ' << tau <<
' ' << xoffset <<
' '
160 << freq <<
' '<< argument
163 std::cerr << retVal << std::endl;
166 else if (_type==
MIXED) {
167 double retVal = (expG-asymm)/2.0;
170 <<
"Warning in AnalyticConvolution: negative probablity"
174 << xsigma <<
' ' << tau <<
' ' << xoffset <<
' '
175 << freq <<
' ' << argument
178 std::cerr << retVal << std::endl;
186 <<
"Unknown sign parity. State is not allowed"
194 double AnalyticConvolution::erfc(
double x)
const {
198 z = (x < 0) ? -x : x;
200 ans = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
201 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
202 t*(-0.82215223+t*0.17087277 ))) ))) )));
203 if ( x < 0 ) ans = 2.0 - ans;
207 double AnalyticConvolution::pow(
double x,
int n)
const {
209 for(
int i=0;i<
n;i++){
215 std::complex<double> AnalyticConvolution::nwwerf(std::complex<double> z)
const {
216 std::complex<double> zh,r[38],s,t,v;
219 const double hf = z1/2;
220 const double z10 = 10;
221 const double c1 = 74/z10;
222 const double c2 = 83/z10;
223 const double c3 = z10/32;
224 const double c4 = 16/z10;
225 const double c = 1.12837916709551257;
226 const double p = pow(2.0*c4,33);
230 double xa=(x >= 0) ? x : -x;
231 double ya=(y >= 0) ? y : -y;
232 if(ya < c1 && xa < c2){
233 zh = std::complex<double>(ya+c4,xa);
234 r[37]= std::complex<double>(0,0);
236 for(
int n = 36; n>0; n--){
237 t=zh+
double(n)*std::conj(r[n+1]);
241 s=std::complex<double>(0,0);
243 for(
int k=33; k>0; k--){
250 zh=std::complex<double>(ya,xa);
251 r[1]=std::complex<double>(0,0);
253 for(
int n=9;n>0;n--){
254 t=zh+
double(n)*std::conj(r[1]);
259 if(ya == 0) v= std::complex<double>(exp(-xa*xa),v.imag());
261 v=2.0*std::exp(std::complex<double>(-xa,-ya)*std::complex<double>(xa,ya))-v;
262 if(x > 0) v=std::conj(v);
265 if(x < 0) v=std::conj(v);
double norm(const HepGenMatrix &m)
AnalyticConvolution(Type=SMEARED_EXP)
virtual double operator()(double argument) const
#define FUNCTION_OBJECT_IMP(classname)
virtual ~AnalyticConvolution()
virtual double getValue() const